COD_20230301_HOST_1_QC_analysis

Minsik Kim

2023-03-14

Loading packages

#===============================================================================
#BTC.LineZero.Header.1.1.0
#===============================================================================
#R Markdown environment setup and reporting utility.
#===============================================================================
#RLB.Dependencies:
#   knitr, magrittr, pacman, rio, rmarkdown, rmdformats, tibble, yaml
#===============================================================================
#Input for document parameters, libraries, file paths, and options.
#=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
knitr::opts_chunk$set(message=FALSE, warning = FALSE)

path_working <- "/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/5_Scripts/MGK/Host_depletion_git/"
path_library <- "/Library/Frameworks/R.framework/Resources/library"
str_libraries <- c(
    "readxl", "phyloseq", "tidyverse", "pacman", "yaml"
)

path_working <- "/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/5_Scripts/MGK/Host_depletion_git"
path_library <- "/Library/Frameworks/R.framework/Resources/library"
str_libraries <- c("readxl", "phyloseq", "tidyverse", "pacman", "yaml", "ggplot2", "vegan", "microbiome", "ggpubr", "viridis", "decontam", "gridExtra", "ggpubr", "lme4", "lmerTest", "writexl", "harrietr", "Maaslin2", "ggtext", "ggpmisc", "gridExtra", "gamm4", "reshape2", "kableExtra", "knitr")
        
YAML_header <-
'---
title: "Host-DNA depletion 1: data wrangling"
author: "Minsik Kim"
date: "2032.03.10"
output:
    rmdformats::downcute:
        downcute_theme: "chaos"
        code_folding: hide
        fig_width: 6
        fig_height: 6
---'
seed <- "20230314"

#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#Loads libraries, file paths, and other document options.
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
FUN.LineZero.Boot <- function() {
    .libPaths(path_library)

    require(pacman)
    pacman::p_load(c("knitr", "rmarkdown", "rmdformats", "yaml"))

    knitr::opts_knit$set(root.dir = path_working)

    str_libraries |> unique() |> sort() -> str_libraries
    pacman::p_load(char = str_libraries)

    set.seed(seed)
}
<<<<<<< HEAD
FUN.LineZero.Boot()
## 
## The downloaded binary packages are in
##  /var/folders/z2/kwm76s1n3jj2sjznjyj9thqc0000gn/T//Rtmp7e3TsQ/downloaded_packages
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
=======
FUN.LineZero.Boot()
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
#Outputs R environment report.
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
FUN.LineZero.Report <- function() {
    cat("Line Zero Environment:\n\n")
    paste("R:", pacman::p_version(), "\n") |> cat()
    cat("Libraries:\n")
    for (str_libraries in str_libraries) {
        paste(
            "    ", str_libraries, ": ", pacman::p_version(package = str_libraries),
            "\n", sep = ""
        ) |> cat()
    }
    paste("\nOperating System:", pacman::p_detectOS(), "\n") |> cat()
    paste("    Library Path:", path_library, "\n") |> cat()
    paste("    Working Path:", path_working, "\n") |> cat()
    paste("Seed:", seed, "\n\n") |> cat()
    cat("YAML Header:\n")
    cat(YAML_header)
}
FUN.LineZero.Report()
## Line Zero Environment:
## 
## R: 4.2.2 
## Libraries:
##     readxl: 1.4.2
<<<<<<< HEAD
##     phyloseq: 1.42.0
=======
##     phyloseq: 1.40.0
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
##     tidyverse: 2.0.0
##     pacman: 0.5.1
##     yaml: 2.3.7
##     ggplot2: 3.4.1
##     vegan: 2.6.4
<<<<<<< HEAD
##     microbiome: 1.20.0
##     ggpubr: 0.6.0
##     viridis: 0.6.2
##     decontam: 1.18.0
=======
##     microbiome: 1.18.0
##     ggpubr: 0.6.0
##     viridis: 0.6.2
##     decontam: 1.16.0
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
##     gridExtra: 2.3
##     ggpubr: 0.6.0
##     lme4: 1.1.31
##     lmerTest: 3.1.3
##     writexl: 1.4.2
##     harrietr: 0.2.3
<<<<<<< HEAD
##     Maaslin2: 1.12.0
=======
##     Maaslin2: 1.10.0
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
##     ggtext: 0.1.2
##     ggpmisc: 0.5.2
##     gridExtra: 2.3
##     gamm4: 0.2.6
##     reshape2: 1.4.4
##     kableExtra: 1.3.4
##     knitr: 1.42
## 
## Operating System: Darwin 
##     Library Path: /Library/Frameworks/R.framework/Resources/library 
##     Working Path: /Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/5_Scripts/MGK/Host_depletion_git 
## Seed: 20230314 
## 
## YAML Header:
## ---
## title: "Host-DNA depletion 1: data wrangling"
## author: "Minsik Kim"
## date: "2032.03.10"
## output:
##     rmdformats::downcute:
##         downcute_theme: "chaos"
##         code_folding: hide
##         fig_width: 6
##         fig_height: 6
## ---

#Script description#

*1. Loading data 1.1. phyloseq obejct 1.2. qPCR data (controls)

  1. QC Q1. How many samples failed sequencing Q2. How were changes in read stats and host DNA proportion? Q3. How were the extraction controls Q4. Prevalence / abundance filtering - red flag

  2. Exploratory analysis *

data input required:

Meta data

  • qPCR - bacteria

  • qPCR - human

  • qPCR host %

  • Raw reads

  • final reads

  • sequencing host %

  • library prep failure status

  • Raw reads

  • subject_id

  • treatment

  • sample_type

  • subject_id

Sequencing result

  • samples

  • controls

Analysis

  1. calculation of alpha-diversity indices

  2. qPCR, by smaple type and treatment A. Total DNA B. Host DNA C. Bacterial DNA D. Host %

  3. Sequencing results (treatment, raw reads, host reads, % host, final reads) –> both stratified and nonstratified library prep % vs sample_type + treatment + sample_type * treatment + subject_id (binary) log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution) Host DNA % vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution) –> both stratified and nonstratified

  4. LM of taxa alpha diversity (richness) –> both stratified and nonstratified

  5. Permanova (Taxa dist ~ log10(final reads) + sample_type + treatment + sample_type * treatment + subject_id) –> both stratified and nonstratified

  6. DA analysis for taxa, by sample type and treatment –> both stratified and nonstratified –> look at level groups as well - phylum, etc.

  7. Decontam - stratified by treatment –> both stratified and nonstratified

  8. LM of taxa alpha diversity (Shannon) –> both stratified and nonstratified

  9. LM of taxa alpha diversity (InvSimpson) –> both stratified and nonstratified

  10. LM of taxa alpha diversity (BPI) –> both stratified and nonstratified

  11. LM of function alpha diversity (richness) –> both stratified and nonstratified

  12. LM of function alpha diversity (Shannon) –> both stratified and nonstratified

  13. LM of function alpha diversity (InvSimpson) –> both stratified and nonstratified

  14. LM of function alpha diversity (BPI) –> both stratified and nonstratified

  15. permanova of function alpha diversity –> both stratified and nonstratified

#Loading data#

# Loading files -----------------------------------------------------------
#loading tidy phyloseq object
phyloseq <- read_rds("/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/4_Data/2_Tidy/Phyloseq/PHY_20221129_MGK_host_tidy_tax.rds")
                
        ####Need to be removed after adding sequenicing data
                data_qPCR <- read_csv("/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/7_Manuscripts/2022_MGK_Host_Depletion/Tables/DAT_20230122_MGK_host_control_qPCR.csv")
                #qPCR data of all the samples sent out sequencing
                data_qPCR <- subset(data_qPCR, data_qPCR$baylor_other_id %in% c(sample_names(phyloseq$phyloseq_count)) | data_qPCR$baylor_other_id %in% c(read_excel("/Users/minsikkim/Dropbox (Partners HealthCare)/Project_Baylor/3_Documentation/Communications/2023-01-24_baylor_shipping_sicas2_nasal_host_depleted/CMMR_MetadataCapture_20230124_LaiP-PQ00430_SICAS2_NS.xlsx", skip = 27) %>% data.frame %>% .$`Optional..............secondary.ID`))
                
                data_qPCR <- subset(data_qPCR, data_qPCR$sample_type %in% c("BAL", "nasal_swab", "Sputum", "neg_depletion", "pos_depletion"))
                data_qPCR$sample_type <- factor(data_qPCR$sample_type, levels = c("BAL", "nasal_swab", "Sputum", "pos_depletion", "neg_depletion"),
                                                labels = c("BAL", "Nasal swab", "Sputum", "P depletion", "N depletion"))
                data_qPCR$treatment <- factor(data_qPCR$treatment, levels = c("control", "lyPMA", "benzonase", "host_zero", "molysis", "qiaamp"),
                                                labels = c("Control", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"))

#sample data loading
sample_data <- sample_data(phyloseq$phyloseq_count)
  1. QC Q1. How many samples failed sequencing Q2. How were changes in read stats and host DNA proportion? Q3. How were the extraction controls Q4. Decontam - were there any contaminants?

Q1. How many samples failed in sequencing?

# Initail QC --------------------------------------------------------------
        #Quesetions - QC
        
#Q0. How many samples failed in sequencing

## figures -----raw data

sample_data %>% 
        subset(., !is.na(.$subject_id)) %>%
        data.frame() %>%
        mutate(host_seq_percent = sequencing_host_prop * 100, 
               .after = sequencing_host_prop,) %>% 
        gather(feature, value, Raw_reads:host_seq_percent) %>%
        group_by(feature, sample_type) %>% 
        subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent")) %>%
        mutate(feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
        ggplot(aes(x = value, fill = treatment)) +
                geom_histogram(bins = 97) +
                theme(legend.position = "top") +
                guides(fill=guide_legend(title="Treatment", nrow = 1)) +
                facet_grid(sample_type~feature, scales = "free") +
                ggtitle("log10 transfromed histrogram")

## figures -----log10

sample_data %>% 
        subset(., !is.na(.$subject_id)) %>%
        data.frame() %>%
        mutate(host_seq_percent = sqrt(sequencing_host_prop), 
               .after = sequencing_host_prop,) %>% 
        gather(feature, value, Raw_reads:host_seq_percent) %>%
        group_by(feature, sample_type) %>% 
        subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent")) %>%
        mutate(feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
        ggplot(aes(x = log10(value), fill = treatment)) +
                geom_histogram(bins = 97) +
                theme(legend.position = "top") +
                guides(fill=guide_legend(title="Treatment", nrow = 1)) +
                facet_grid(sample_type~feature, scales = "free") +
                ggtitle("log10 for reads, sqrt for host % transfromed histrogram")

# no somaple showed 0 reads

ggarrange(ggplot(sample_data %>% subset(., !is.na(.$subject_id)) %>% data.frame(), aes(x = Final_reads, fill = treatment)) +
                geom_histogram(bins = 97) +
                facet_wrap(~sample_type) +
                theme_classic(base_family = "serif") +
                ggtitle("Histogram of final reads by sample type and treatment") +
                scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), name = "Treatment"),
             ggplot(sample_data %>% subset(., !is.na(.$subject_id)) %>% data.frame(), aes(x = log10(Final_reads), fill = treatment)) +
                geom_histogram(bins = 97) +
                facet_wrap(~sample_type) +
                theme_classic(base_family = "serif") +
                ggtitle("Histogram of log10(final reads) by sample type and treatment") +
                scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), name = "Treatment"),
          common.legend = T, ncol = 1)

hist((log10((phyloseq$phyloseq_count %>% otu_table %>% colSums()) + 1)),100, main = "Histogram of total reads (sum of OTU table)") # 2 samples showed 0 total reads (sum of otu_table)

#how were the samples failed in library prep?
sample_data %>% data.frame %>% mutate(total_read = phyloseq$phyloseq_count %>% otu_table %>% colSums()) %>% rownames_to_column(var = "baylor_other_id") %>%
        ggplot(aes(x = reorder(baylor_other_id, -total_read),
                               y = log10(total_read + 1),
                               col = lib_failed)) +
                geom_point() +
                theme_classic(base_family = "serif") +
                theme(axis.title.y = element_markdown(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 4)) +
                ylab("log<sub>10</sub>(Sum of OTU table reads)") +
                xlab("Sample ID") +
        guides(col=guide_legend(title="Library failed")) +
        ggtitle("Sum of OTU reads by library failure status")

sample_data %>% data.frame %>% filter(phyloseq$phyloseq_count %>% otu_table %>% colSums() == 0) # two BAL sampels showed 0 total reads
sample_data(phyloseq$phyloseq_count) %>% data.frame() %>% subset(., .$lib_failed)

Results:

2.1.1 Modeling final read should be conducted with both linear and logistic regression. Host % need no transformation, and should be tested for both linear and binomial methods

2.1.2 13 samples failed in library prep

2.1.3. Two BAL sampels showed 0 total reads

2.1.4. Sequencing fail ≠ library prep fail

Comments from Baylor:

Q: What was the lab’s criteria for deciding which samples failed library prep.? There were 13 samples that you pointed as failed but their sequencing result actually looks pretty good (ie similar to samples that didn’t fail library prep)

A: To determine whether a library attempt “passed or failed” the lab looks at the picogreen concentrations and a library fragment size distribution trace. The trace files are an output from either the Fragment Analyzer or TapeStation (a copy of the trace files for PQ00331 is attached). If a sample has a background level pico concentration and no discernable fragment concentration on the trace (i.e. a flat trace line) it is considered failed library. If the sample is below the level of detection of our standard library QC methods, it is considered failure. It’s still possible that there is some small amounts of library in those samples that were successfully sequenced, but often those samples do not generate a meaningful amount of sequence data.

QC2 table

#sequencing result by sample type and control (1/0)
options(dplyr.summarise.inform = FALSE)

sample_data %>% data.frame() %>% 
        group_by(sample_type, treated) %>% 
        summarise(N = n(),
                  `Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2),nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
        ) %>%
        rename(`Sample type` = sample_type, Treated = treated) %>%
        data.frame(check.names = F) %>% mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")
Sample type Treated N Raw reads
(median [IQR])
[reads x 107] 		Host reads
(median [IQR])
[reads x 107] 		Host reads proportion
(median [IQR])
[%] 		Final reads
(median [IQR])
[reads x 107]
BAL 0 5 15.73 [6.35, 15.92] 12.92 [5.21, 12.94] 82.14 [82.12, 82.18] 0.03 [0.03, 0.04]
BAL 1 25 6.17 [4.57, 17.43] 4.65 [2.78, 12.80] 75.59 [69.38, 78.44] 0.17 [0.10, 0.37]
Nasal swab 0 10 13.09 [7.73, 16.93] 10.05 [6.11, 13.04] 77.34 [76.89, 79.45] 0.48 [0.10, 0.87]
Nasal swab 1 25 4.08 [0.99, 6.40] 0.81 [0.26, 1.36] 27.46 [13.51, 64.58] 0.97 [0.17, 3.42]
Sputum 0 5 8.59 [8.25, 9.27] 6.87 [6.69, 7.50] 80.61 [79.92, 80.89] 0.06 [0.06, 0.09]
Sputum 1 25 12.23 [10.34, 13.73] 7.71 [3.76, 8.82] 58.64 [39.71, 74.93] 1.16 [0.47, 4.19]
Neg. ext. 0 1 0.02 [0.02, 0.02] 0.00 [0.00, 0.00] 3.13 [3.13, 3.13] 0.02 [0.02, 0.02]
Pos. ext. 0 1 11.87 [11.87, 11.87] 0.00 [0.00, 0.00] 0.02 [0.02, 0.02] 9.96 [9.96, 9.96] 
sample_data %>% data.frame() %>% 
        #dplyr::filter(sample_type %in% c("Sputum", "nasal_swab", "BAL")) %>% 
        group_by (sample_type, treatment) %>%
        summarise(N = n(),
              `Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
              `Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
              `Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
              `Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
        ) %>% data.frame(check.names = F) %>% 
        arrange(sample_type, treatment) %>%
        rename(`Sample type` = sample_type, Treatment = treatment) %>%
        mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")
Sample type Treatment N Raw reads
(median [IQR])
[reads x 107] 		Host reads
(median [IQR])
[reads x 107] 		Host reads proportion
(median [IQR])
[%] 		Final reads
(median [IQR])
[reads x 107]
BAL Control 5 15.73 [6.35, 15.92] 12.92 [5.21, 12.94] 82.14 [82.12, 82.18] 0.03 [0.03, 0.04]
BAL lyPMA 5 5.72 [3.59, 13.41] 4.65 [2.79, 10.90] 80.73 [77.85, 81.18] 0.06 [0.04, 0.10]
BAL Benzonase 5 18.59 [16.20, 23.63] 14.77 [12.80, 18.16] 79.02 [77.89, 79.45] 0.17 [0.16, 0.22]
BAL Host zero 5 4.57 [2.32, 4.71] 2.69 [1.61, 2.93] 66.98 [57.17, 68.95] 0.24 [0.13, 0.82]
BAL Molysis 5 4.76 [3.57, 4.86] 2.78 [1.39, 3.61] 74.39 [74.29, 75.59] 0.29 [0.13, 1.56]
BAL QIAamp 5 17.19 [15.35, 17.43] 11.87 [10.79, 12.22] 74.88 [71.08, 77.31] 0.26 [0.10, 1.02]
Nasal swab Control 10 13.09 [7.73, 16.93] 10.05 [6.11, 13.04] 77.34 [76.89, 79.45] 0.48 [0.10, 0.87]
Nasal swab lyPMA 5 0.98 [0.85, 1.24] 0.63 [0.28, 0.88] 71.37 [28.73, 74.68] 0.07 [0.06, 0.08]
Nasal swab Benzonase 5 5.75 [4.95, 6.57] 3.66 [1.29, 5.05] 64.58 [63.71, 76.83] 0.28 [0.26, 1.04]
Nasal swab Host zero 5 2.83 [1.42, 6.42] 0.49 [0.03, 0.81] 7.67 [2.33, 25.73] 2.43 [0.97, 5.03]
Nasal swab Molysis 5 0.99 [0.63, 4.08] 0.42 [0.06, 0.64] 41.04 [4.31, 64.24] 0.32 [0.17, 2.53]
Nasal swab QIAamp 5 6.40 [6.40, 6.80] 0.86 [0.86, 1.17] 17.24 [13.51, 19.57] 4.63 [4.50, 4.67]
Sputum Control 5 8.59 [8.25, 9.27] 6.87 [6.69, 7.50] 80.61 [79.92, 80.89] 0.06 [0.06, 0.09]
Sputum lyPMA 5 10.98 [5.22, 12.78] 8.82 [3.76, 10.44] 76.33 [72.02, 80.29] 0.25 [0.15, 0.44]
Sputum Benzonase 5 10.76 [10.34, 10.82] 7.81 [7.75, 8.24] 74.93 [72.13, 75.33] 0.47 [0.45, 0.59]
Sputum Host zero 5 13.14 [7.64, 13.95] 4.39 [3.80, 7.71] 49.71 [31.45, 55.49] 2.91 [2.36, 3.67]
Sputum Molysis 5 12.59 [10.84, 13.73] 2.98 [1.83, 4.28] 27.49 [14.62, 28.19] 6.11 [5.56, 8.37]
Sputum QIAamp 5 12.35 [12.23, 12.85] 9.08 [8.41, 9.27] 71.76 [56.69, 73.50] 1.16 [1.13, 3.89]
Neg. ext. Control 1 0.02 [0.02, 0.02] 0.00 [0.00, 0.00] 3.13 [3.13, 3.13] 0.02 [0.02, 0.02]
Pos. ext. Control 1 11.87 [11.87, 11.87] 0.00 [0.00, 0.00] 0.02 [0.02, 0.02] 9.96 [9.96, 9.96]

QC2 figure

# Summary figures - facet and z-score -------------------------------------

sample_data %>% 
        subset(., !is.na(.$subject_id)) %>%
        data.frame() %>%
        gather(feature, value, Raw_reads:sequencing_host_prop) %>%
        group_by(feature, sample_type) %>% 
        subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop")) %>%
        mutate(z_score = scale(value),
               feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
        ggplot(aes(x = treatment, y = z_score, fill = treatment)) +
                geom_boxplot(lwd = 0.2) +
                guides(fill=guide_legend(title="Treatment", nrow = 1)) +
                facet_grid(sample_type~feature) +
                xlab("Treatment") +
                ylab("Z score") +
                theme_classic(base_family = "serif", base_size = 14) +
                guides( x =  guide_axis(angle = 90)) + 
                theme(legend.position = "top") +
                scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), name = "Treatment") #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6

Results:

2.2.1. There were no differences between raw reads, however, final reads increased after some treatment, and host DNA proportion decreased

2.2.2. The changes in final reads were varied by treatment and sample type

comment: Should I normalize z-score to mean value of control?

Q3. How was the extraction controls, both pos and neg?

#Loading theoretical mock community

zymo_mock <- read_excel("/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_sicas2/data_raw/DAR_20210929_zymo_mock_data.xlsx") %>%
        data.frame(row.names = T) %>% rename(mock_theoretical = Mock) %>% mutate(mock_theoretical = mock_theoretical/100) %>%
        merge_phyloseq(otu_table(., taxa_are_rows = T), tax_table(phyloseq$phyloseq_count))

phyloseq_control_rel <- rbind(c("mock_theoretical", "mock")) %>% data.frame() %>%
        column_to_rownames(var = "X1") %>% rename(sample_type = X2) %>% #making sample_data of mock community
        merge_phyloseq(sample_data(.), zymo_mock) %>%        
        merge_phyloseq(., subset_samples(phyloseq$phyloseq_rel, sample_type == "Pos. ext." | sample_type == "Neg. ext." )) #adding data of controls


#Species richness of each control groups
cat("Species richness\n")
## Species richness
rowSums(t(otu_table(phyloseq_control_rel)) != 0)
## mock_theoretical     20220606_Pos     20220606_Neg 
##               10               41                3
sample_data(phyloseq_control_rel)$sample_type <- factor(sample_data(phyloseq_control_rel)$sample_type,
                                                        levels = c("mock", "2", "1"),
                                                        labels = c("Mock, 10 spp.", "Pos., 41 spp.", "Neg., 3 spp."))

#Manipulating phyloseq - only top 10 

tax_table(phyloseq_control_rel) %>%
  cbind(species10 = "[Others]") %>%
  {top10species <- head(taxa_sums(phyloseq_control_rel), 10) %>% names()
   .[top10species, "species10"] <- as.character(.[top10species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species10") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Species richness of control data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ sample_type, scales= "free_x", nrow=1)

#there could be opportunistic pathogens...

Results

2.3.1. Negative control showed minimal number of possible contaminants

2.3.2. Positive control contained various contaminants

Results:

2.2.1. There were no differences between raw reads, however, final reads increased after some treatment, and host DNA proportion decreased

2.2.2. The changes in final reads were varied by treatment and sample type

Q4. Prevalence filtering - red flag

#Calculation of sample prevalence, standard deviation, median abundance across all samples for all bugs and making into a table.
#
#•  In initial analysis we will not perform prevalence or abundance filtering (though we may consider this for secondary differential abundance analyses to manage p (features) > n (sample size) problem and issues with multiple hypothesis correction)
taxa_qc <- data.frame("species" =  otu_table(phyloseq$phyloseq_rel) %>% t() %>% colnames(),
        "prevalence" = ifelse(phyloseq$phyloseq_count %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
        "mean_rel_abd" = phyloseq$phyloseq_rel %>% otu_table() %>% t() %>% colMeans(na.rm = T) #mean relativ abundacne 
)


hist(log10(taxa_qc$prevalence), xlab = "log10(Taxa prevalence)", main = "Histogram of prevalence of taxa")

hist(log10(taxa_qc$mean_rel_abd), xlab = "log10(Mean relative abundance)", main = "Histogram of mean relative abundance")

• In initial analysis we will not perform prevalence or abundance filtering (though we may consider this for secondary differentialabundance analyses to manage p (features) > n (sample size) problem and issues with multiple hypothesis correction) • Calculation of sample prevalence, standard deviation, median abundance across all samples for all bugs and making into a table.

Although we don’t consider the prevalence of abundance at this time, we can consider their red-flags after running the DA analysis

    1. red flag
red_flag_taxa <- data.frame(species = taxa_qc$species,
                          red_flag_prev_abd = ifelse(taxa_qc$prevalence < otu_table(phyloseq$phyloseq_rel) %>% t %>% rownames() %>% length * 0.05 & taxa_qc$mean_rel_abd < quantile(taxa_qc$mean_rel_abd, 0.75), 1, 0)
)

#Analysis#

A0. calculation of alpha-diversity indices

alpha_diversity <- function(data) {
        otu_table <- otu_table(data) %>% .[colSums(.) !=0]
        S.obs <- rowSums(t(otu_table) != 0)
        sample_data <- sample_data(data)
        data_evenness <- vegan::diversity(t(otu_table)) / log(vegan::specnumber(t(otu_table))) # calculate evenness index using vegan package
        data_shannon <- vegan::diversity(t(otu_table), index = "shannon") # calculate Shannon index using vegan package
        data_hill <- exp(data_shannon)                           # calculate Hills index
        
        data_dominance <- microbiome::dominance(otu_table, index = "all", rank = 1, aggregate = TRUE) # dominance (Berger-Parker index), etc.
        data_invsimpson <- vegan::diversity(t(otu_table), index = "invsimpson")                          # calculate Shannon index using vegan package
        alpha_diversity <- cbind(S.obs, data_shannon, data_hill, data_invsimpson, data_evenness,data_dominance) # combine all indices in one data table
        sample_data <- merge(data.frame(sample_data), alpha_diversity, by = 0, all = T) %>% column_to_rownames(var = "Row.names")
}

sample_data(phyloseq$phyloseq_rel) <- sample_data(alpha_diversity(phyloseq$phyloseq_count))
sample_data(phyloseq$phyloseq_count) <- sample_data(alpha_diversity(phyloseq$phyloseq_count)) 
sample_data(phyloseq$phyloseq_path_rpkm) <- sample_data(alpha_diversity(phyloseq$phyloseq_path_rpkm))

A1. qPCR, by smaple type and treatment

    A. Total DNA
    B. Host DNA
    C. Bacterial DNA
    D. Host %
#2A: Change in total DNA (qPCR)


f2a <- ggplot(data_qPCR, aes(x = sample_type, y = log10(DNA_host_nondil + DNA_bac_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        ylab("log<sub>10</sub>(qPCR total DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif") + 
        labs(tag = "A") +
        scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))


#2B: Change in human DNA (qPCR)
f2b <- ggplot(data_qPCR, aes(x = sample_type, y = log10(DNA_host_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) +
        ylab("log<sub>10</sub>(qPCR host DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif")+ 
        labs(tag = "B") +
        scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#2C: Change in 16S DNA (qPCR)
f2c <- ggplot(data_qPCR, aes(x = sample_type, y = log10(DNA_bac_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) +
        ylab("log<sub>10</sub>(qPCR bacterial DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif")+ 
        labs(tag = "C") +
        scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))

#2D. Change in % host (qPCR)
f2d <- ggplot(data_qPCR, aes(x = sample_type, y = log10(host_proportion * 100))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) +
        ylab("log<sub>10</sub>(host DNA) (%)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif") + 
        labs(tag = "D") +
        scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))

#output for markdown
ggarrange(f2a, f2b, f2c, f2d, common.legend = T , align = "hv")

A2. Sequencing results (treatment, raw reads, host reads, % host, final reads)

--> both stratified and nonstratified
 library prep % vs sample_type + treatment + sample_type * treatment + subject_id (binary)
log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution)
Host DNA % vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution)
--> both stratified and nonstratified
cat("Binomial - library failed? \n\n")
## Binomial - library failed?
glmer(lib_failed ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error t value
(Intercept) 0.0000000 0.1167011 0.0000000
sample_typeNasal swab -0.0204626 0.1470919 -0.1391145
sample_typeSputum 0.0000000 0.1650403 0.0000000
treatmentlyPMA 0.2000000 0.1619675 1.2348160
treatmentBenzonase 0.0000000 0.1619675 0.0000000
treatmentHost zero 0.2000000 0.1619675 1.2348160
treatmentMolysis 0.2000000 0.1619675 1.2348160
treatmentQIAamp 0.0000000 0.1619675 0.0000000
sample_typeNasal swab:treatmentlyPMA 0.5978181 0.2143680 2.7887474
sample_typeSputum:treatmentlyPMA -0.2000000 0.2290566 -0.8731468
sample_typeNasal swab:treatmentBenzonase -0.0021819 0.2143680 -0.0101781
sample_typeSputum:treatmentBenzonase 0.0000000 0.2290566 0.0000000
sample_typeNasal swab:treatmentHost zero 0.1978181 0.2143680 0.9227971
sample_typeSputum:treatmentHost zero -0.2000000 0.2290566 -0.8731468
sample_typeNasal swab:treatmentMolysis 0.6021819 0.2143680 2.8091035
sample_typeSputum:treatmentMolysis -0.2000000 0.2290566 -0.8731468
sample_typeNasal swab:treatmentQIAamp 0.0021819 0.2143680 0.0101781
sample_typeSputum:treatmentQIAamp 0.0000000 0.2290566 0.0000000 
cat("Binomial - Final reads? \n\n")
## Binomial - Final reads?
glmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error t value
(Intercept) 5.4933063 0.2095210 26.2184086
sample_typeNasal swab 1.1246781 0.2789651 4.0316080
sample_typeSputum 0.3332900 0.2963074 1.1248117
treatmentlyPMA 0.3521385 0.2733402 1.2882794
treatmentBenzonase 0.8109049 0.2733402 2.9666509
treatmentHost zero 0.9501904 0.2733402 3.4762193
treatmentMolysis 1.0358039 0.2733402 3.7894313
treatmentQIAamp 1.0480985 0.2733402 3.8344107
sample_typeNasal swab:treatmentlyPMA -0.8774152 0.3621899 -2.4225278
sample_typeSputum:treatmentlyPMA 0.1879364 0.3865614 0.4861747
sample_typeNasal swab:treatmentBenzonase -0.6887274 0.3621899 -1.9015642
sample_typeSputum:treatmentBenzonase 0.0349655 0.3865614 0.0904526
sample_typeNasal swab:treatmentHost zero -0.1157121 0.3621899 -0.3194791
sample_typeSputum:treatmentHost zero 0.7231403 0.3865614 1.8706997
sample_typeNasal swab:treatmentMolysis -0.8411412 0.3621899 -2.3223760
sample_typeSputum:treatmentMolysis 0.9543008 0.3865614 2.4686913
sample_typeNasal swab:treatmentQIAamp 0.0236242 0.3621899 0.0652260
sample_typeSputum:treatmentQIAamp 0.3676901 0.3865614 0.9511817 
cat("Linear - Final reads? \n\n")
## Linear - Final reads?
lmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.4933063 0.2095210 62.79820 26.2184086 0.0000000
sample_typeNasal swab 1.1246781 0.2789651 37.26869 4.0316080 0.0002633
sample_typeSputum 0.3332900 0.2963074 62.79820 1.1248117 0.2649498
treatmentlyPMA 0.3521385 0.2733402 67.65016 1.2882794 0.2020376
treatmentBenzonase 0.8109049 0.2733402 67.65016 2.9666509 0.0041580
treatmentHost zero 0.9501904 0.2733402 67.65016 3.4762193 0.0008929
treatmentMolysis 1.0358039 0.2733402 67.65016 3.7894313 0.0003238
treatmentQIAamp 1.0480985 0.2733402 67.65016 3.8344107 0.0002788
sample_typeNasal swab:treatmentlyPMA -0.8774152 0.3621899 68.18953 -2.4225278 0.0180749
sample_typeSputum:treatmentlyPMA 0.1879364 0.3865614 67.65016 0.4861747 0.6284145
sample_typeNasal swab:treatmentBenzonase -0.6887274 0.3621899 68.18953 -1.9015642 0.0614539
sample_typeSputum:treatmentBenzonase 0.0349655 0.3865614 67.65016 0.0904526 0.9281949
sample_typeNasal swab:treatmentHost zero -0.1157121 0.3621899 68.18953 -0.3194791 0.7503399
sample_typeSputum:treatmentHost zero 0.7231403 0.3865614 67.65016 1.8706997 0.0657128
sample_typeNasal swab:treatmentMolysis -0.8411412 0.3621899 68.18953 -2.3223760 0.0232052
sample_typeSputum:treatmentMolysis 0.9543008 0.3865614 67.65016 2.4686913 0.0160928
sample_typeNasal swab:treatmentQIAamp 0.0236242 0.3621899 68.18953 0.0652260 0.9481849
sample_typeSputum:treatmentQIAamp 0.3676901 0.3865614 67.65016 0.9511817 0.3448982 
cat("Binomial - Host DNA proportion? \n\n")
## Binomial - Host DNA proportion?
glmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error t value
(Intercept) 0.8200713 0.0646272 12.6892671
sample_typeNasal swab -0.0378253 0.0799149 -0.4733196
sample_typeSputum -0.0293956 0.0913966 -0.3216265
treatmentlyPMA -0.0377605 0.0908962 -0.4154243
treatmentBenzonase -0.0324363 0.0908962 -0.3568498
treatmentHost zero -0.2000025 0.0908962 -2.2003393
treatmentMolysis -0.1577230 0.0908962 -1.7351990
treatmentQIAamp -0.0928027 0.0908962 -1.0209737
sample_typeNasal swab:treatmentlyPMA -0.1931415 0.1202628 -1.6059959
sample_typeSputum:treatmentlyPMA 0.0106590 0.1285467 0.0829196
sample_typeNasal swab:treatmentBenzonase -0.1301969 0.1202628 -1.0826036
sample_typeSputum:treatmentBenzonase -0.0432386 0.1285467 -0.3363648
sample_typeNasal swab:treatmentHost zero -0.3924910 0.1202628 -3.2636117
sample_typeSputum:treatmentHost zero -0.1532698 0.1285467 -1.1923278
sample_typeNasal swab:treatmentMolysis -0.2637315 0.1202628 -2.1929608
sample_typeSputum:treatmentMolysis -0.3863225 0.1285467 -3.0053090
sample_typeNasal swab:treatmentQIAamp -0.5240678 0.1202628 -4.3576894
sample_typeSputum:treatmentQIAamp -0.0370057 0.1285467 -0.2878775 
cat("Linear - Host DNA proportion? \n\n")
## Linear - Host DNA proportion?
lmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.8200713 0.0646272 75.67093 12.6892671 0.0000000
sample_typeNasal swab -0.0378253 0.0799149 41.00023 -0.4733196 0.6384957
sample_typeSputum -0.0293956 0.0913966 75.67093 -0.3216265 0.7486224
treatmentlyPMA -0.0377605 0.0908962 68.23467 -0.4154243 0.6791348
treatmentBenzonase -0.0324363 0.0908962 68.23467 -0.3568498 0.7223055
treatmentHost zero -0.2000025 0.0908962 68.23467 -2.2003393 0.0311715
treatmentMolysis -0.1577230 0.0908962 68.23467 -1.7351990 0.0872202
treatmentQIAamp -0.0928027 0.0908962 68.23467 -1.0209737 0.3108732
sample_typeNasal swab:treatmentlyPMA -0.1931415 0.1202628 68.80454 -1.6059959 0.1128550
sample_typeSputum:treatmentlyPMA 0.0106590 0.1285467 68.23467 0.0829196 0.9341582
sample_typeNasal swab:treatmentBenzonase -0.1301969 0.1202628 68.80454 -1.0826036 0.2827639
sample_typeSputum:treatmentBenzonase -0.0432386 0.1285467 68.23467 -0.3363648 0.7376281
sample_typeNasal swab:treatmentHost zero -0.3924910 0.1202628 68.80454 -3.2636117 0.0017148
sample_typeSputum:treatmentHost zero -0.1532698 0.1285467 68.23467 -1.1923278 0.2372628
sample_typeNasal swab:treatmentMolysis -0.2637315 0.1202628 68.80454 -2.1929608 0.0316937
sample_typeSputum:treatmentMolysis -0.3863225 0.1285467 68.23467 -3.0053090 0.0037094
sample_typeNasal swab:treatmentQIAamp -0.5240678 0.1202628 68.80454 -4.3576894 0.0000449
sample_typeSputum:treatmentQIAamp -0.0370057 0.1285467 68.23467 -0.2878775 0.7743131

A3. LM of taxa alpha diversity

--> both stratified and nonstratified
sample_data <- sample_data(phyloseq$phyloseq_count) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))

cat("Species richness of all samples\n\n")
## Species richness of all samples
cat("S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -62.286176 18.946348 73.91996 -3.2875031 0.0015486
sample_typeNasal swab -9.894010 13.143390 19.76996 -0.7527746 0.4604558
sample_typeSputum 5.066489 11.833705 23.07241 0.4281405 0.6725190
treatmentlyPMA -3.032607 7.992068 64.29757 -0.3794521 0.7056025
treatmentBenzonase -4.719045 8.051585 65.30664 -0.5861014 0.5598277
treatmentHost zero -2.064590 8.208256 65.44300 -0.2515260 0.8021956
treatmentMolysis 10.862491 8.314322 65.52729 1.3064795 0.1959566
treatmentQIAamp -2.491588 8.330143 65.53938 -0.2991050 0.7658063
log10(Final_reads) 12.532134 3.113656 67.56866 4.0248934 0.0001465
sample_typeNasal swab:treatmentlyPMA 4.202591 10.419319 64.56848 0.4033461 0.6880263
sample_typeSputum:treatmentlyPMA 34.664316 10.599221 64.21630 3.2704588 0.0017287
sample_typeNasal swab:treatmentBenzonase 1.775047 10.042806 64.92798 0.1767481 0.8602564
sample_typeSputum:treatmentBenzonase 62.518484 10.357810 64.45145 6.0358784 0.0000001
sample_typeNasal swab:treatmentHost zero 4.793943 9.787379 64.68743 0.4898086 0.6259264
sample_typeSputum:treatmentHost zero 92.094185 10.559614 64.43191 8.7213589 0.0000000
sample_typeNasal swab:treatmentMolysis -7.289176 10.178107 65.15977 -0.7161623 0.4764500
sample_typeSputum:treatmentMolysis 87.197251 10.723044 64.48528 8.1317631 0.0000000
sample_typeNasal swab:treatmentQIAamp -3.726532 9.774308 64.66718 -0.3812579 0.7042616
sample_typeSputum:treatmentQIAamp 70.548734 10.400892 64.40670 6.7829501 0.0000000 
#Species richness  - stratified 
cat("Species richness at NS\n\n")
## Species richness at NS
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal swab")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -26.0233860 9.340261 28 -2.7861520 0.0094654
treatmentlyPMA -2.1153025 2.324810 28 -0.9098818 0.3706512
treatmentBenzonase -1.7637014 2.206148 28 -0.7994486 0.4307605
treatmentHost zero 8.8224892 2.495655 28 3.5351399 0.0014386
treatmentMolysis 4.5783098 2.217883 28 2.0642703 0.0483678
treatmentQIAamp 0.8360832 2.678401 28 0.3121575 0.7572336
log10(Final_reads) 5.6349921 1.419898 28 3.9685903 0.0004571 
cat("Species richness at BAL\n\n")
## Species richness at BAL
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -100.395549 25.504894 20.38095 -3.9363249 0.0007925
treatmentlyPMA -5.839094 7.201885 17.39942 -0.8107731 0.4284473
treatmentBenzonase -10.667780 7.752673 18.42945 -1.3760131 0.1853100
treatmentHost zero -8.986746 8.091470 18.61512 -1.1106444 0.2808621
treatmentMolysis 3.342010 8.316699 18.72132 0.4018433 0.6923498
treatmentQIAamp -10.097992 8.350027 18.73603 -1.2093364 0.2415732
log10(Final_reads) 19.520813 4.535759 19.97726 4.3037587 0.0003466 
cat("Species richness at sputum\n\n")
## Species richness at sputum
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -166.73637 94.19800 20.78908 -1.770063 0.0913826
treatmentlyPMA 21.48046 13.48708 19.68002 1.592670 0.1271695
treatmentBenzonase 41.90046 17.06102 19.99448 2.455918 0.0233254
treatmentHost zero 58.57768 28.77613 20.32431 2.035634 0.0550337
treatmentMolysis 60.65374 33.57099 20.37080 1.806731 0.0855997
treatmentQIAamp 41.44599 24.96239 20.26633 1.660337 0.1122410
log10(Final_reads) 31.32813 16.05007 20.49824 1.951899 0.0647464 
#Shannon
cat("Shannon index\n\n")
## Shannon index
cat("Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_shannon ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.8948018 0.8689906 73.16148 2.1804629 0.0324397
sample_typeNasal swab 0.3612177 0.4227092 52.85894 0.8545299 0.3966679
sample_typeSputum 0.5229775 0.4502631 46.00549 1.1614932 0.2514323
treatmentlyPMA 0.1432302 0.3457153 59.27054 0.4143009 0.6801493
treatmentBenzonase 0.3245806 0.3510837 61.59790 0.9245107 0.3588282
treatmentHost zero 0.1622396 0.3591607 62.11606 0.4517187 0.6530450
treatmentMolysis 0.9220443 0.3646163 62.43914 2.5288070 0.0139853
treatmentQIAamp 0.3416438 0.3654292 62.48560 0.9349109 0.3534333
log10(Final_reads) -0.2037567 0.1481509 70.28159 -1.3753323 0.1733965
sample_typeNasal swab:treatmentlyPMA -0.3349724 0.4621995 62.11790 -0.7247356 0.4713354
sample_typeSputum:treatmentlyPMA 0.9305639 0.4573666 58.90852 2.0346129 0.0463989
sample_typeNasal swab:treatmentBenzonase -0.1716183 0.4431174 62.20574 -0.3872975 0.6998579
sample_typeSputum:treatmentBenzonase 1.1940106 0.4466493 59.19803 2.6732618 0.0096921
sample_typeNasal swab:treatmentHost zero 0.1836222 0.4309190 61.71186 0.4261176 0.6715048
sample_typeSputum:treatmentHost zero 1.4844081 0.4573881 59.56239 3.2454017 0.0019263
sample_typeNasal swab:treatmentMolysis -0.6161157 0.4488821 62.42148 -1.3725557 0.1748022
sample_typeSputum:treatmentMolysis 0.9186726 0.4659713 59.93997 1.9715216 0.0532853
sample_typeNasal swab:treatmentQIAamp -0.2466370 0.4305635 61.76297 -0.5728238 0.5688436
sample_typeSputum:treatmentQIAamp 1.0739999 0.4489918 59.21646 2.3920256 0.0199485 
#Simpson
cat("Inverse Simpson index\n\n")
## Inverse Simpson index
cat("InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_invsimpson ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.0016985 4.1346231 73.75264 1.2097109 0.2302531
sample_typeNasal swab -0.1347032 1.9113437 60.99582 -0.0704757 0.9440456
sample_typeSputum 0.3376540 2.0144643 55.69325 0.1676148 0.8674935
treatmentlyPMA -0.0446283 1.6994301 60.26753 -0.0262607 0.9791362
treatmentBenzonase -0.0159936 1.7179035 63.20493 -0.0093099 0.9926012
treatmentHost zero -0.5969192 1.7555995 63.84335 -0.3400088 0.7349661
treatmentMolysis 2.8388719 1.7811019 64.24055 1.5938852 0.1158705
treatmentQIAamp -0.2830161 1.7849044 64.29760 -0.1585609 0.8745118
log10(Final_reads) -0.4205707 0.7102118 72.89587 -0.5921764 0.5555648
sample_typeNasal swab:treatmentlyPMA -0.2501880 2.2597896 63.56198 -0.1107130 0.9121928
sample_typeSputum:treatmentlyPMA 3.6784976 2.2498385 59.81795 1.6350052 0.1073001
sample_typeNasal swab:treatmentBenzonase 0.5135753 2.1660273 63.71628 0.2371047 0.8133362
sample_typeSputum:treatmentBenzonase 8.4907847 2.1958416 60.21080 3.8667565 0.0002731
sample_typeNasal swab:treatmentHost zero 1.1519699 2.1086029 63.04119 0.5463190 0.5867750
sample_typeSputum:treatmentHost zero 7.4853477 2.2470527 60.66550 3.3311847 0.0014760
sample_typeNasal swab:treatmentMolysis -2.0541758 2.1931699 64.02401 -0.9366241 0.3524716
sample_typeSputum:treatmentMolysis 4.7103881 2.2875456 61.13466 2.0591450 0.0437470
sample_typeNasal swab:treatmentQIAamp 0.4812870 2.1067028 63.07313 0.2284551 0.8200315
sample_typeSputum:treatmentQIAamp 5.7791599 2.2072771 60.23515 2.6182303 0.0111632 
#BPI
cat("Beger Parker index\n\n")
## Beger Parker index
cat("BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(dbp ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.0429030 0.3089340 73.34596 -0.1388743 0.8899301
sample_typeNasal swab -0.2254845 0.1470717 55.46609 -1.5331600 0.1309210
sample_typeSputum -0.1615566 0.1559725 48.85446 -1.0358016 0.3053947
treatmentlyPMA -0.0992353 0.1244133 59.03815 -0.7976258 0.4282842
treatmentBenzonase -0.2090856 0.1261407 61.68022 -1.6575576 0.1024850
treatmentHost zero -0.1071198 0.1289953 62.26501 -0.8304161 0.4094764
treatmentMolysis -0.3789439 0.1309245 62.62965 -2.8943694 0.0052227
treatmentQIAamp -0.2002048 0.1312121 62.68208 -1.5258106 0.1320873
log10(Final_reads) 0.1455197 0.0528497 71.25889 2.7534627 0.0074764
sample_typeNasal swab:treatmentlyPMA 0.1074786 0.1660074 62.17039 0.6474324 0.5197358
sample_typeSputum:treatmentlyPMA -0.1999551 0.1646342 58.63042 -1.2145416 0.2294101
sample_typeNasal swab:treatmentBenzonase 0.0783117 0.1591430 62.28479 0.4920838 0.6243901
sample_typeSputum:treatmentBenzonase -0.2242934 0.1607441 58.96684 -1.3953444 0.1681450
sample_typeNasal swab:treatmentHost zero -0.0843810 0.1548174 61.70491 -0.5450359 0.5876935
sample_typeSputum:treatmentHost zero -0.4349677 0.1645678 59.37801 -2.6430911 0.0104899
sample_typeNasal swab:treatmentMolysis 0.1755526 0.1611877 62.54326 1.0891189 0.2802803
sample_typeSputum:treatmentMolysis -0.2575968 0.1676125 59.80380 -1.5368593 0.1296021
sample_typeNasal swab:treatmentQIAamp 0.1096753 0.1546846 61.75206 0.7090255 0.4809770
sample_typeSputum:treatmentQIAamp -0.2873906 0.1615851 58.98802 -1.7785712 0.0804633

A4. Taxa beta diversity

Permanova (Taxa dist ~ log10(final reads) + sample_type + treatment + sample_type * treatment + subject_id) –> both stratified and nonstratified

phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_rel, S.obs != 0 & sample_type %in% c("BAL", "Nasal swab", "Sputum"))

bray_perm_uni <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + treatment + subject_id,
                            data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)

bray_perm <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp + subject_id,
                            data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)

bray_perm_strata <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp,
                            data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F),
                            strata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F) %>% .$subject_id, permutations = 10000)

bray_perm_inter <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type * treatment + log10(Final_reads) + subject_id,
                                  data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

bray_perm_ns <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
                               data = subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab") %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

bray_perm_bal  <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "BAL"), method="bray") ~  lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
                                 data = subset_samples(phyloseq_rel_nz, sample_type == "BAL") %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

bray_perm_spt <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
                                data = subset_samples(phyloseq_rel_nz, sample_type == "Sputum") %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

cat("\nUnivariate analysis\n")
## 
## Univariate analysis
bray_perm_uni %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 13.429 0.371 56.251 0.000
log10(Final reads) 1 1.111 0.031 9.310 0.000
Treatment 5 1.213 0.033 2.033 0.001
Subject 10 11.639 0.321 9.751 0.000
Residual 74 8.833 0.244 NA NA
Total 92 36.225 1.000 NA NA 
cat("\nDetailed treatment\n")
## 
## Detailed treatment
bray_perm %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 13.429 0.371 56.251 0.000
log10(Final reads) 1 1.111 0.031 9.310 0.000
lyPMA 1 0.140 0.004 1.172 0.287
Benzonase 1 0.092 0.003 0.774 0.607
Host zero 1 0.157 0.004 1.312 0.218
Molysis 1 0.200 0.006 1.676 0.103
QIAamp 1 0.624 0.017 5.229 0.000
Subject 10 11.639 0.321 9.751 0.000
Residual 74 8.833 0.244 NA NA
Total 92 36.225 1.000 NA NA 
cat("\n Strata -detailed treatment\n")
## 
##  Strata -detailed treatment
bray_perm_strata %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 13.429 0.371 27.550 0.000
log10(Final reads) 1 1.111 0.031 4.560 0.000
lyPMA 1 0.140 0.004 0.574 0.530
Benzonase 1 0.092 0.003 0.379 0.665
Host zero 1 0.157 0.004 0.643 0.406
Molysis 1 0.200 0.006 0.821 0.167
QIAamp 1 0.624 0.017 2.561 0.003
Residual 84 20.472 0.565 NA NA
Total 92 36.225 1.000 NA NA 
cat("\n Interaction term \n")
## 
##  Interaction term
bray_perm_inter %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "sample_type:treatment" ~ 'Sample type X treatment',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 13.429 0.371 65.605 0
Treatment 5 1.184 0.033 2.314 0
log10(Final reads) 1 1.140 0.031 11.142 0
Subject 10 11.639 0.321 11.372 0
Sample type X treatment 10 2.283 0.063 2.231 0
Residual 64 6.550 0.181 NA NA
Total 92 36.225 1.000 NA NA 
cat("\n Stratified - nasal swab \n")
## 
##  Stratified - nasal swab
bray_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
lyPMA 1 0.730 0.120 5.122 0.004
Benzonase 1 0.191 0.031 1.342 0.262
Host zero 1 0.171 0.028 1.201 0.306
Molysis 1 0.137 0.022 0.961 0.405
QIAamp 1 0.254 0.042 1.780 0.148
log10(Final reads) 1 0.428 0.070 3.007 0.036
Subject id 2 0.488 0.080 1.714 0.119
Residual 26 3.704 0.607 NA NA
Total 34 6.103 1.000 NA NA 
cat("\n Stratified - BAL \n")
## 
##  Stratified - BAL
bray_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
lyPMA 1 0.100 0.010 1.144 0.308
Benzonase 1 0.025 0.003 0.288 0.934
Host zero 1 0.086 0.009 0.990 0.418
Molysis 1 0.085 0.009 0.969 0.428
QIAamp 1 0.229 0.024 2.625 0.035
log10(Final reads) 1 1.482 0.152 16.966 0.000
Subject id 4 6.241 0.641 17.864 0.000
Residual 17 1.485 0.153 NA NA
Total 27 9.734 1.000 NA NA 
cat("\n Stratified - sputum \n")
## 
##  Stratified - sputum
bray_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
lyPMA 1 0.139 0.020 2.565 0.046
Benzonase 1 0.037 0.005 0.688 0.635
Host zero 1 0.171 0.025 3.150 0.018
Molysis 1 0.436 0.063 8.046 0.000
QIAamp 1 0.953 0.137 17.588 0.000
log10(Final reads) 1 0.172 0.025 3.174 0.017
Subject id 4 4.022 0.578 18.555 0.000
Residual 19 1.030 0.148 NA NA
Total 29 6.960 1.000 NA NA

A5. DA analysis for taxa, by sample type and treatment

--> both stratified and nonstratified
--> look at level groups as well - phylum, etc.
<<<<<<< HEAD 
#DA analysis - MaAslin

sample_data(phyloseq_rel_nz)$log10.Final_reads <- log10(sample_data(phyloseq_rel_nz)$Final_reads)
sample_data(phyloseq_rel_nz)$sampletype_treatment <- paste(sample_data(phyloseq_rel_nz)$sample_type, sample_data(phyloseq_rel_nz)$treatment, sep = ":")

#Running MaAslin for all sample without decontam
#for taxa differentially abundant by host depletion method, look to see which ones overlap with potential contaminant taxa

# Maaslin - # # y ~ log(final reads) + sample_type + treatment  -----------

#all samples
maaslin_all = Maaslin2(input_data = otu_table(phyloseq_rel_nz) %>% t %>% data.frame(), 
                input_metadata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), 
                output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output", 
                fixed_effects = c("sample_type", "log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"), 
                transform = "LOG", #default
                normalization = "TSS", # default
                random_effects = c("subject_id"), 
                reference = c("sample_type,BAL"), 
                plot_heatmap = F,
                plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:26 INFO::Writing function arguments to log file
## 2023-03-22 00:45:26 INFO::Verifying options selected are valid
## 2023-03-22 00:45:26 INFO::Determining format of input files
## 2023-03-22 00:45:26 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:26 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:26 INFO::Formula for fixed effects: expr ~  sample_type + log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:26 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:26 INFO::Total samples in data: 93
## 2023-03-22 00:45:26 INFO::Min samples required with min abundance for a feature not to be filtered: 9.300000
## 2023-03-22 00:45:26 INFO::Total filtered features: 183
## 2023-03-22 00:45:26 INFO::Filtered feature names from abundance and prevalence filtering: Actinomyces_sp_oral_taxon_414, Actinomyces_sp_oral_taxon_448, Actinomyces_sp_oral_taxon_897, Parascardovia_denticolens, Corynebacterium_afermentans, Corynebacterium_aurimucosum, Corynebacterium_kroppenstedtii, Corynebacterium_matruchotii, Corynebacterium_tuberculostearicum, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Tropheryma_whipplei, Cutibacterium_avidum, Atopobium_deltae, Adlercreutzia_equolifaciens, Rubrobacter_radiotolerans, Patulibacter_medicamentivorans, Prevotella_oris, Capnocytophaga_leadbetteri, Brochothrix_campestris, Staphylococcus_capitis, Enterococcus_faecalis, Lactobacillus_rhamnosus, Streptococcus_sinensis, Eubacterium_nodatum, Lachnospiraceae_bacterium_oral_taxon_096, Oribacterium_parvum, Shuttleworthia_satelles, Selenomonas_noxia, Selenomonas_sputigena, Anaerococcus_nagyae, Anaerococcus_octavius, Peptoniphilus_lacrimalis, Leptotrichia_wadei, Achromobacter_denitrificans, Achromobacter_ruhlandii, Achromobacter_xylosoxidans, Lautropia_mirabilis, Leptothrix_ochracea, Cardiobacterium_hominis, Escherichia_coli, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Fretibacterium_fastidiosum, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_tropicalis, Malassezia_globosa, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_radingae, Actinomyces_turicensis, Varibaculum_cambriense, Bifidobacterium_breve, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Nakamurella_silvestris, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_asaccharolytica, Porphyromonas_canoris, Porphyromonas_catoniae, Porphyromonas_uenonis, Alloprevotella_rava, Alloprevotella_tannerae, Prevotella_buccae, Prevotella_denticola, Prevotella_intermedia, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pleuritidis, Prevotella_scopos, Prevotella_sp_F0091, Prevotella_sp_oral_taxon_306, Prevotella_veroralis, Tannerella_forsythia, Tannerella_sp_oral_taxon_808, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Litorilinea_aerophila, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_fermentum, Lactobacillus_gasseri, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Sharpea_azabuensis, Streptococcus_massiliensis, Streptococcus_sp_SK643, Streptococcus_sp_oral_taxon_056, Mageeibacillus_indolicus, Oribacterium_asaccharolyticum, Peptostreptococcus_sp_MV1, Leptotrichia_sp_oral_taxon_212, Desulfobulbus_oralis, Moraxella_catarrhalis, Stenotrophomonas_pavanii, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Staphylococcus_haemolyticus, Enterococcus_avium, Streptococcus_downei, Streptococcus_pyogenes, Streptococcus_salivarius_CAG_79, Streptococcus_sobrinus, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_thermophilus, Streptococcus_viridans, Eubacterium_saphenum, Eubacterium_rectale, Filifactor_alocis, Eggerthia_catenaformis, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_invisus, Dialister_micraerophilus, Negativicoccus_succinicivorans, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_215, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Haemophilus_parainfluenzae, Candida_orthopsilosis, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_concisus, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Serratia_marcescens, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:26 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:26 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:26 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:26 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:26 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:26 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 1, Actinobaculum_sp_oral_taxon_183
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 2, Actinomyces_georgiae
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 3, Actinomyces_graevenitzii
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 4, Actinomyces_hongkongensis
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 5, Actinomyces_johnsonii
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 6, Actinomyces_massiliensis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 7, Actinomyces_meyeri
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 8, Actinomyces_naeslundii
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 9, Actinomyces_odontolyticus
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 10, Actinomyces_oris
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 11, Actinomyces_sp_HMSC035G02
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 12, Actinomyces_sp_HPA0247
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 13, Actinomyces_sp_ICM47
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 14, Actinomyces_sp_S6_Spd3
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 15, Actinomyces_sp_oral_taxon_170
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 16, Actinomyces_sp_oral_taxon_180
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 17, Actinomyces_sp_oral_taxon_181
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 18, Actinomyces_viscosus
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 19, Aeriscardovia_aeriphila
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 20, Alloscardovia_omnicolens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 21, Bifidobacterium_dentium
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 22, Bifidobacterium_longum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 23, Scardovia_wiggsiae
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 24, Corynebacterium_accolens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 25, Corynebacterium_atypicum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 26, Corynebacterium_durum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 27, Corynebacterium_pseudodiphtheriticum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 28, Corynebacterium_pseudogenitalium
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 29, Rothia_aeria
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 30, Rothia_dentocariosa
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 31, Rothia_mucilaginosa
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 32, Cutibacterium_acnes
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 33, Cutibacterium_granulosum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 34, Propionibacterium_humerusii
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 35, Propionibacterium_namnetense
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 36, Propionibacterium_acidifaciens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 37, Pseudopropionibacterium_propionicum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 38, Atopobium_parvulum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 39, Atopobium_rimae
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 40, Olsenella_profusa
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 41, Olsenella_scatoligenes
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 42, Olsenella_uli
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 43, Collinsella_aerofaciens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 44, Collinsella_intestinalis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 45, Collinsella_stercoris
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 46, Enorma_massiliensis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 47, X.Collinsella._massiliensis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 48, Asaccharobacter_celatus
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 49, Denitrobacterium_detoxificans
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 50, Enterorhabdus_caecimuris
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 51, Slackia_exigua
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 52, Slackia_isoflavoniconvertens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 53, Thermoleophilum_album
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 54, Porphyromonas_endodontalis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 55, Porphyromonas_somerae
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 56, Prevotella_histicola
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 57, Prevotella_melaninogenica
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 58, Kouleothrix_aurantiaca
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 59, Hydrogenibacillus_schlegelii
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 60, Gemella_asaccharolytica
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 61, Gemella_bergeri
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 62, Gemella_haemolysans
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 63, Gemella_morbillorum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 64, Gemella_sanguinis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 65, Staphylococcus_argenteus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 66, Staphylococcus_aureus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 67, Staphylococcus_epidermidis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 68, Staphylococcus_schweitzeri
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 69, Abiotrophia_defectiva
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 70, Abiotrophia_sp_HMSC24B09
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 71, Dolosigranulum_pigrum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 72, Granulicatella_adiacens
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 73, Granulicatella_elegans
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 74, Streptococcus_anginosus_group
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 75, Streptococcus_australis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 76, Streptococcus_cristatus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 77, Streptococcus_gordonii
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 78, Streptococcus_infantis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 79, Streptococcus_milleri
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 80, Streptococcus_mitis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 81, Streptococcus_oralis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 82, Streptococcus_parasanguinis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 83, Streptococcus_peroris
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 84, Streptococcus_pneumoniae
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 85, Streptococcus_pseudopneumoniae
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 86, Streptococcus_salivarius
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 87, Streptococcus_sanguinis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 88, Streptococcus_sp_A12
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 89, Streptococcus_sp_F0442
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 90, Streptococcus_sp_HMSC034E03
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 91, Streptococcus_sp_HMSC067H01
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 92, Streptococcus_sp_HMSC071D03
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 93, Streptococcus_sp_HPH0090
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 94, Streptococcus_sp_M334
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 95, Streptococcus_vestibularis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 96, Eubacterium_brachy
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 97, Eubacterium_infirmum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 98, Eubacterium_sulci
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 99, Mogibacterium_diversum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 100, Mogibacterium_pumilum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 101, Mogibacterium_timidum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 102, Lachnoanaerobaculum_saburreum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 103, Oribacterium_sinus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 104, Oribacterium_sp_oral_taxon_078
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 105, Stomatobaculum_longum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 106, Peptostreptococcus_stomatis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 107, Bulleidia_extructa
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 108, Solobacterium_moorei
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 109, Limnochorda_pilosa
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 110, Veillonella_atypica
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 111, Veillonella_dispar
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 112, Veillonella_infantium
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 113, Veillonella_parvula
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 114, Veillonella_sp_T11011_6
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 115, Finegoldia_magna
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 116, Parvimonas_micra
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 117, Parvimonas_sp_oral_taxon_110
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 118, Parvimonas_sp_oral_taxon_393
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 119, Fusobacterium_nucleatum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 120, Gemmata_obscuriglobus
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 121, Paludisphaera_borealis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 122, Cupriavidus_sp
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 123, Sutterella_parvirubra
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 124, Anaerobiospirillum_thomasii
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 125, Cardiobacterium_valvarum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 126, Alkalilimnicola_ehrlichii
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 127, Thiohalorhabdus_denitrificans
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 128, Pseudomonas_aeruginosa_group
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 129, Acholeplasma_oculi
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 130, Candida_albicans
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 131, Candida_dubliniensis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 132, Candida_parapsilosis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 133, Malassezia_restricta
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 134, Prevotella_jejuni
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 135, Prevotella_pallens
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 136, Prevotella_salivae
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 137, Tannerella_sp_oral_taxon_HOT_286
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 138, Capnocytophaga_gingivalis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 139, Streptococcus_mutans
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 140, Megasphaera_micronuciformis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 141, Neisseria_flavescens
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 142, Neisseria_subflava
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 143, Stenotrophomonas_maltophilia
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 144, Stenotrophomonas_rhizophila
## 2023-03-22 00:45:29 INFO::Counting total values for each feature
## 2023-03-22 00:45:29 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:29 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:29 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:29 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:29 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:29 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:29 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:29 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
#Checking number of bugs for the main text of the manuscript 
cat("Number of differentially abundant bugs by each metadata")
## Number of differentially abundant bugs by each metadata
maaslin_all$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                 8                 6                50                14 
##           molysis            qiaamp       sample_type 
##                 9                 2                66
cat("Decreased bugs by each metadata")
## Decreased bugs by each metadata
maaslin_all$results %>% subset(., .$qval < 0.1 & .$coef < 0) %>% .$metadata %>% table()
## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                 1                 1                 1                 2 
##           molysis       sample_type 
##                 1                 1
cat("Incrased")
## Incrased
maaslin_all$results %>% subset(., .$qval < 0.1 & .$coef > 0) %>% .$metadata %>% table()
## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                 7                 5                49                12 
##           molysis            qiaamp       sample_type 
##                 8                 2                65
#NS
# # y ~ log(final reads) + sample_type + treatment 

fit_data_ns = Maaslin2(input_data = otu_table(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab")) %>% t %>% data.frame(),
                       input_metadata = sample_data(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab")) %>% data.frame(), 
                       output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output", 
                       fixed_effects = c("log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"), 
                       transform = "LOG", #default
                       normalization = "TSS", # default
                       random_effects = c("subject_id"), 
                       plot_heatmap = F,
                       plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:29 INFO::Writing function arguments to log file
## 2023-03-22 00:45:29 INFO::Verifying options selected are valid
## 2023-03-22 00:45:29 INFO::Determining format of input files
## 2023-03-22 00:45:29 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:29 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:29 INFO::Formula for fixed effects: expr ~  log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:29 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:29 INFO::Total samples in data: 35
## 2023-03-22 00:45:29 INFO::Min samples required with min abundance for a feature not to be filtered: 3.500000
## 2023-03-22 00:45:29 INFO::Total filtered features: 295
## 2023-03-22 00:45:29 INFO::Filtered feature names from abundance and prevalence filtering: Actinobaculum_sp_oral_taxon_183, Actinomyces_georgiae, Actinomyces_graevenitzii, Actinomyces_hongkongensis, Actinomyces_johnsonii, Actinomyces_massiliensis, Actinomyces_meyeri, Actinomyces_naeslundii, Actinomyces_odontolyticus, Actinomyces_oris, Actinomyces_sp_HMSC035G02, Actinomyces_sp_HPA0247, Actinomyces_sp_ICM47, Actinomyces_sp_S6_Spd3, Actinomyces_sp_oral_taxon_170, Actinomyces_sp_oral_taxon_180, Actinomyces_sp_oral_taxon_181, Actinomyces_sp_oral_taxon_414, Actinomyces_sp_oral_taxon_448, Actinomyces_sp_oral_taxon_897, Actinomyces_viscosus, Alloscardovia_omnicolens, Bifidobacterium_dentium, Bifidobacterium_longum, Parascardovia_denticolens, Scardovia_wiggsiae, Corynebacterium_afermentans, Corynebacterium_durum, Corynebacterium_kroppenstedtii, Corynebacterium_matruchotii, Corynebacterium_tuberculostearicum, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Rothia_aeria, Rothia_dentocariosa, Tropheryma_whipplei, Propionibacterium_acidifaciens, Pseudopropionibacterium_propionicum, Atopobium_deltae, Atopobium_parvulum, Atopobium_rimae, Olsenella_profusa, Olsenella_scatoligenes, Olsenella_uli, Collinsella_stercoris, Enorma_massiliensis, Adlercreutzia_equolifaciens, Asaccharobacter_celatus, Denitrobacterium_detoxificans, Enterorhabdus_caecimuris, Slackia_exigua, Rubrobacter_radiotolerans, Patulibacter_medicamentivorans, Porphyromonas_endodontalis, Porphyromonas_somerae, Prevotella_histicola, Prevotella_melaninogenica, Prevotella_oris, Capnocytophaga_leadbetteri, Kouleothrix_aurantiaca, Hydrogenibacillus_schlegelii, Gemella_bergeri, Gemella_morbillorum, Brochothrix_campestris, Staphylococcus_schweitzeri, Abiotrophia_defectiva, Abiotrophia_sp_HMSC24B09, Granulicatella_adiacens, Granulicatella_elegans, Enterococcus_faecalis, Lactobacillus_rhamnosus, Streptococcus_anginosus_group, Streptococcus_australis, Streptococcus_cristatus, Streptococcus_gordonii, Streptococcus_infantis, Streptococcus_milleri, Streptococcus_mitis, Streptococcus_parasanguinis, Streptococcus_peroris, Streptococcus_pneumoniae, Streptococcus_pseudopneumoniae, Streptococcus_salivarius, Streptococcus_sanguinis, Streptococcus_sinensis, Streptococcus_sp_A12, Streptococcus_sp_F0442, Streptococcus_sp_HMSC034E03, Streptococcus_sp_HMSC067H01, Streptococcus_sp_HMSC071D03, Streptococcus_sp_HPH0090, Streptococcus_sp_M334, Streptococcus_vestibularis, Eubacterium_brachy, Eubacterium_infirmum, Eubacterium_nodatum, Eubacterium_sulci, Mogibacterium_diversum, Mogibacterium_pumilum, Mogibacterium_timidum, Lachnoanaerobaculum_saburreum, Lachnospiraceae_bacterium_oral_taxon_096, Oribacterium_parvum, Oribacterium_sinus, Oribacterium_sp_oral_taxon_078, Shuttleworthia_satelles, Stomatobaculum_longum, Peptostreptococcus_stomatis, Bulleidia_extructa, Solobacterium_moorei, Selenomonas_noxia, Selenomonas_sputigena, Veillonella_atypica, Veillonella_dispar, Veillonella_infantium, Veillonella_parvula, Veillonella_sp_T11011_6, Anaerococcus_octavius, Parvimonas_micra, Parvimonas_sp_oral_taxon_110, Parvimonas_sp_oral_taxon_393, Peptoniphilus_lacrimalis, Fusobacterium_nucleatum, Leptotrichia_wadei, Gemmata_obscuriglobus, Paludisphaera_borealis, Achromobacter_denitrificans, Achromobacter_ruhlandii, Achromobacter_xylosoxidans, Lautropia_mirabilis, Leptothrix_ochracea, Anaerobiospirillum_thomasii, Cardiobacterium_hominis, Cardiobacterium_valvarum, Alkalilimnicola_ehrlichii, Escherichia_coli, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Thiohalorhabdus_denitrificans, Pseudomonas_aeruginosa_group, Fretibacterium_fastidiosum, Acholeplasma_oculi, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_albicans, Candida_dubliniensis, Candida_parapsilosis, Candida_tropicalis, Malassezia_globosa, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_radingae, Actinomyces_turicensis, Varibaculum_cambriense, Bifidobacterium_breve, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Nakamurella_silvestris, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_asaccharolytica, Porphyromonas_canoris, Porphyromonas_catoniae, Porphyromonas_uenonis, Alloprevotella_rava, Alloprevotella_tannerae, Prevotella_buccae, Prevotella_denticola, Prevotella_intermedia, Prevotella_jejuni, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pallens, Prevotella_pleuritidis, Prevotella_salivae, Prevotella_scopos, Prevotella_sp_F0091, Prevotella_sp_oral_taxon_306, Prevotella_veroralis, Tannerella_forsythia, Tannerella_sp_oral_taxon_808, Tannerella_sp_oral_taxon_HOT_286, Capnocytophaga_gingivalis, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Litorilinea_aerophila, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_fermentum, Lactobacillus_gasseri, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Sharpea_azabuensis, Streptococcus_massiliensis, Streptococcus_mutans, Streptococcus_sp_SK643, Streptococcus_sp_oral_taxon_056, Mageeibacillus_indolicus, Oribacterium_asaccharolyticum, Peptostreptococcus_sp_MV1, Megasphaera_micronuciformis, Leptotrichia_sp_oral_taxon_212, Neisseria_flavescens, Neisseria_subflava, Desulfobulbus_oralis, Moraxella_catarrhalis, Stenotrophomonas_maltophilia, Stenotrophomonas_pavanii, Stenotrophomonas_rhizophila, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Staphylococcus_haemolyticus, Enterococcus_avium, Streptococcus_downei, Streptococcus_pyogenes, Streptococcus_salivarius_CAG_79, Streptococcus_sobrinus, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_thermophilus, Streptococcus_viridans, Eubacterium_saphenum, Eubacterium_rectale, Filifactor_alocis, Eggerthia_catenaformis, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_invisus, Dialister_micraerophilus, Negativicoccus_succinicivorans, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_215, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Haemophilus_parainfluenzae, Candida_orthopsilosis, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_concisus, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Serratia_marcescens, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:29 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:29 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:29 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:29 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:29 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:29 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 1, Aeriscardovia_aeriphila
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 2, Corynebacterium_accolens
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 3, Corynebacterium_atypicum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 4, Corynebacterium_aurimucosum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 5, Corynebacterium_pseudodiphtheriticum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 6, Corynebacterium_pseudogenitalium
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 7, Rothia_mucilaginosa
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 8, Cutibacterium_acnes
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 9, Cutibacterium_avidum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 10, Cutibacterium_granulosum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 11, Propionibacterium_humerusii
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 12, Propionibacterium_namnetense
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 13, Collinsella_aerofaciens
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 14, Collinsella_intestinalis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 15, X.Collinsella._massiliensis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 16, Slackia_isoflavoniconvertens
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 17, Thermoleophilum_album
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 18, Gemella_asaccharolytica
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 19, Gemella_haemolysans
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 20, Gemella_sanguinis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 21, Staphylococcus_argenteus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 22, Staphylococcus_aureus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 23, Staphylococcus_capitis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 24, Staphylococcus_epidermidis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 25, Dolosigranulum_pigrum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 26, Streptococcus_oralis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 27, Limnochorda_pilosa
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 28, Anaerococcus_nagyae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 29, Finegoldia_magna
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 30, Cupriavidus_sp
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 31, Sutterella_parvirubra
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 32, Malassezia_restricta
## 2023-03-22 00:45:30 INFO::Counting total values for each feature
## 2023-03-22 00:45:30 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:30 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:30 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:30 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:30 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:30 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:30 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:30 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
#Checking number of bugs for the main text of the manuscript 
cat("Number of differentially abundant bugs by each metadata in NS")
## Number of differentially abundant bugs by each metadata in NS
fit_data_ns$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                 2                 6                 9                 4 
##           molysis            qiaamp 
##                 4                 1
#bal
# # y ~ log(final reads) + sample_type + treatment 

fit_data_bal = Maaslin2(input_data = otu_table(subset_samples(phyloseq_rel_nz, sample_type == "BAL")) %>% t %>% data.frame(),
                       input_metadata = sample_data(subset_samples(phyloseq_rel_nz, sample_type == "BAL")) %>% data.frame(), 
                       output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output", 
                       fixed_effects = c("log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"), 
                       transform = "LOG", #default
                       normalization = "TSS", # default
                       random_effects = c("subject_id"), 
                       plot_heatmap = F,
                       plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:30 INFO::Writing function arguments to log file
## 2023-03-22 00:45:30 INFO::Verifying options selected are valid
## 2023-03-22 00:45:30 INFO::Determining format of input files
## 2023-03-22 00:45:30 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:30 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:30 INFO::Formula for fixed effects: expr ~  log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:30 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:30 INFO::Total samples in data: 28
## 2023-03-22 00:45:30 INFO::Min samples required with min abundance for a feature not to be filtered: 2.800000
## 2023-03-22 00:45:30 INFO::Total filtered features: 261
## 2023-03-22 00:45:30 INFO::Filtered feature names from abundance and prevalence filtering: Actinobaculum_sp_oral_taxon_183, Actinomyces_georgiae, Actinomyces_graevenitzii, Actinomyces_hongkongensis, Actinomyces_johnsonii, Actinomyces_massiliensis, Actinomyces_meyeri, Actinomyces_oris, Actinomyces_sp_HMSC035G02, Actinomyces_sp_ICM47, Actinomyces_sp_oral_taxon_170, Actinomyces_sp_oral_taxon_180, Actinomyces_sp_oral_taxon_414, Actinomyces_sp_oral_taxon_448, Actinomyces_sp_oral_taxon_897, Actinomyces_viscosus, Aeriscardovia_aeriphila, Alloscardovia_omnicolens, Bifidobacterium_dentium, Bifidobacterium_longum, Parascardovia_denticolens, Corynebacterium_accolens, Corynebacterium_afermentans, Corynebacterium_atypicum, Corynebacterium_aurimucosum, Corynebacterium_durum, Corynebacterium_kroppenstedtii, Corynebacterium_matruchotii, Corynebacterium_pseudodiphtheriticum, Corynebacterium_pseudogenitalium, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Rothia_aeria, Rothia_dentocariosa, Tropheryma_whipplei, Cutibacterium_avidum, Cutibacterium_granulosum, Propionibacterium_humerusii, Propionibacterium_namnetense, Propionibacterium_acidifaciens, Pseudopropionibacterium_propionicum, Olsenella_profusa, Olsenella_uli, Collinsella_aerofaciens, Collinsella_intestinalis, Collinsella_stercoris, Enorma_massiliensis, X.Collinsella._massiliensis, Adlercreutzia_equolifaciens, Asaccharobacter_celatus, Denitrobacterium_detoxificans, Enterorhabdus_caecimuris, Slackia_exigua, Slackia_isoflavoniconvertens, Rubrobacter_radiotolerans, Patulibacter_medicamentivorans, Thermoleophilum_album, Porphyromonas_endodontalis, Porphyromonas_somerae, Prevotella_histicola, Capnocytophaga_leadbetteri, Kouleothrix_aurantiaca, Hydrogenibacillus_schlegelii, Brochothrix_campestris, Staphylococcus_capitis, Abiotrophia_defectiva, Abiotrophia_sp_HMSC24B09, Dolosigranulum_pigrum, Granulicatella_adiacens, Granulicatella_elegans, Lactobacillus_rhamnosus, Streptococcus_australis, Streptococcus_cristatus, Streptococcus_gordonii, Streptococcus_infantis, Streptococcus_milleri, Streptococcus_peroris, Streptococcus_pseudopneumoniae, Streptococcus_salivarius, Streptococcus_sanguinis, Streptococcus_sinensis, Streptococcus_sp_A12, Streptococcus_sp_F0442, Streptococcus_sp_HMSC034E03, Streptococcus_sp_HMSC067H01, Streptococcus_sp_HMSC071D03, Streptococcus_sp_HPH0090, Streptococcus_sp_M334, Streptococcus_vestibularis, Eubacterium_brachy, Eubacterium_nodatum, Eubacterium_sulci, Mogibacterium_diversum, Mogibacterium_pumilum, Mogibacterium_timidum, Lachnoanaerobaculum_saburreum, Lachnospiraceae_bacterium_oral_taxon_096, Oribacterium_parvum, Oribacterium_sinus, Oribacterium_sp_oral_taxon_078, Shuttleworthia_satelles, Stomatobaculum_longum, Solobacterium_moorei, Limnochorda_pilosa, Selenomonas_noxia, Selenomonas_sputigena, Veillonella_atypica, Veillonella_dispar, Veillonella_infantium, Veillonella_sp_T11011_6, Anaerococcus_nagyae, Anaerococcus_octavius, Parvimonas_sp_oral_taxon_110, Parvimonas_sp_oral_taxon_393, Peptoniphilus_lacrimalis, Leptotrichia_wadei, Gemmata_obscuriglobus, Paludisphaera_borealis, Achromobacter_denitrificans, Achromobacter_ruhlandii, Achromobacter_xylosoxidans, Lautropia_mirabilis, Leptothrix_ochracea, Anaerobiospirillum_thomasii, Cardiobacterium_hominis, Alkalilimnicola_ehrlichii, Escherichia_coli, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Thiohalorhabdus_denitrificans, Pseudomonas_aeruginosa_group, Fretibacterium_fastidiosum, Acholeplasma_oculi, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_tropicalis, Malassezia_globosa, Malassezia_restricta, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_turicensis, Bifidobacterium_breve, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Nakamurella_silvestris, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_canoris, Porphyromonas_catoniae, Porphyromonas_uenonis, Prevotella_denticola, Prevotella_intermedia, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pallens, Prevotella_pleuritidis, Prevotella_scopos, Prevotella_sp_F0091, Tannerella_forsythia, Tannerella_sp_oral_taxon_808, Tannerella_sp_oral_taxon_HOT_286, Capnocytophaga_gingivalis, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Litorilinea_aerophila, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_gasseri, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Sharpea_azabuensis, Streptococcus_massiliensis, Streptococcus_mutans, Streptococcus_sp_SK643, Streptococcus_sp_oral_taxon_056, Mageeibacillus_indolicus, Oribacterium_asaccharolyticum, Leptotrichia_sp_oral_taxon_212, Neisseria_subflava, Desulfobulbus_oralis, Moraxella_catarrhalis, Stenotrophomonas_maltophilia, Stenotrophomonas_pavanii, Stenotrophomonas_rhizophila, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Streptococcus_downei, Streptococcus_salivarius_CAG_79, Streptococcus_sobrinus, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_viridans, Eubacterium_saphenum, Eubacterium_rectale, Filifactor_alocis, Eggerthia_catenaformis, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_micraerophilus, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_215, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Haemophilus_parainfluenzae, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:30 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:30 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:30 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:30 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:30 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:30 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 1, Actinomyces_naeslundii
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 2, Actinomyces_odontolyticus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 3, Actinomyces_sp_HPA0247
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 4, Actinomyces_sp_S6_Spd3
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 5, Actinomyces_sp_oral_taxon_181
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 6, Scardovia_wiggsiae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 7, Corynebacterium_tuberculostearicum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 8, Rothia_mucilaginosa
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 9, Cutibacterium_acnes
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 10, Atopobium_deltae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 11, Atopobium_parvulum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 12, Atopobium_rimae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 13, Olsenella_scatoligenes
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 14, Prevotella_melaninogenica
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 15, Prevotella_oris
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 16, Gemella_asaccharolytica
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 17, Gemella_bergeri
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 18, Gemella_haemolysans
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 19, Gemella_morbillorum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 20, Gemella_sanguinis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 21, Staphylococcus_argenteus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 22, Staphylococcus_aureus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 23, Staphylococcus_epidermidis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 24, Staphylococcus_schweitzeri
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 25, Enterococcus_faecalis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 26, Streptococcus_anginosus_group
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 27, Streptococcus_mitis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 28, Streptococcus_oralis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 29, Streptococcus_parasanguinis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 30, Streptococcus_pneumoniae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 31, Eubacterium_infirmum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 32, Peptostreptococcus_stomatis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 33, Bulleidia_extructa
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 34, Veillonella_parvula
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 35, Finegoldia_magna
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 36, Parvimonas_micra
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 37, Fusobacterium_nucleatum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 38, Cupriavidus_sp
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 39, Sutterella_parvirubra
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 40, Cardiobacterium_valvarum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 41, Candida_albicans
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 42, Candida_dubliniensis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 43, Candida_parapsilosis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 44, Actinomyces_radingae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 45, Varibaculum_cambriense
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 46, Porphyromonas_asaccharolytica
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 47, Alloprevotella_rava
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 48, Alloprevotella_tannerae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 49, Prevotella_buccae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 50, Prevotella_jejuni
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 51, Prevotella_salivae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 52, Prevotella_sp_oral_taxon_306
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 53, Prevotella_veroralis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 54, Lactobacillus_fermentum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 55, Peptostreptococcus_sp_MV1
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 56, Megasphaera_micronuciformis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 57, Neisseria_flavescens
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 58, Staphylococcus_haemolyticus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 59, Enterococcus_avium
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 60, Streptococcus_pyogenes
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 61, Streptococcus_thermophilus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 62, Dialister_invisus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 63, Negativicoccus_succinicivorans
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 64, Candida_orthopsilosis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 65, Campylobacter_concisus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 66, Serratia_marcescens
## 2023-03-22 00:45:31 INFO::Counting total values for each feature
## 2023-03-22 00:45:31 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:31 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:31 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:31 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:31 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:31 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:31 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:31 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
cat("Number of differentially abundant bugs by each metadata in BAL")
## Number of differentially abundant bugs by each metadata in BAL
fit_data_bal$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## < table of extent 0 >
#sputum
# # y ~ log(final reads) + sample_type + treatment 

fit_data_spt = Maaslin2(input_data = otu_table(subset_samples(phyloseq_rel_nz, sample_type == "Sputum")) %>% t %>% data.frame(),
                       input_metadata = sample_data(subset_samples(phyloseq_rel_nz, sample_type == "Sputum")) %>% data.frame(), 
                       output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output", 
                       fixed_effects = c("log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"), 
                       transform = "LOG", #default
                       normalization = "TSS", # default
                       random_effects = c("subject_id"),
                       plot_heatmap = F,
                       plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:31 INFO::Writing function arguments to log file
## 2023-03-22 00:45:31 INFO::Verifying options selected are valid
## 2023-03-22 00:45:31 INFO::Determining format of input files
## 2023-03-22 00:45:31 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:31 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:31 INFO::Formula for fixed effects: expr ~  log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:31 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:31 INFO::Total samples in data: 30
## 2023-03-22 00:45:31 INFO::Min samples required with min abundance for a feature not to be filtered: 3.000000
## 2023-03-22 00:45:31 INFO::Total filtered features: 155
## 2023-03-22 00:45:31 INFO::Filtered feature names from abundance and prevalence filtering: Parascardovia_denticolens, Corynebacterium_accolens, Corynebacterium_afermentans, Corynebacterium_aurimucosum, Corynebacterium_kroppenstedtii, Corynebacterium_pseudodiphtheriticum, Corynebacterium_pseudogenitalium, Corynebacterium_tuberculostearicum, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Tropheryma_whipplei, Cutibacterium_avidum, Cutibacterium_granulosum, Propionibacterium_humerusii, Propionibacterium_namnetense, Patulibacter_medicamentivorans, Prevotella_oris, Gemella_asaccharolytica, Staphylococcus_capitis, Dolosigranulum_pigrum, Streptococcus_sinensis, Selenomonas_noxia, Selenomonas_sputigena, Anaerococcus_nagyae, Anaerococcus_octavius, Finegoldia_magna, Lautropia_mirabilis, Leptothrix_ochracea, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_albicans, Candida_dubliniensis, Candida_tropicalis, Malassezia_globosa, Malassezia_restricta, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_radingae, Actinomyces_turicensis, Varibaculum_cambriense, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_asaccharolytica, Porphyromonas_canoris, Porphyromonas_uenonis, Alloprevotella_rava, Alloprevotella_tannerae, Prevotella_buccae, Prevotella_denticola, Prevotella_intermedia, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pleuritidis, Prevotella_scopos, Prevotella_sp_F0091, Prevotella_sp_oral_taxon_306, Prevotella_veroralis, Tannerella_sp_oral_taxon_808, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_fermentum, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Streptococcus_massiliensis, Streptococcus_sp_SK643, Mageeibacillus_indolicus, Peptostreptococcus_sp_MV1, Leptotrichia_sp_oral_taxon_212, Desulfobulbus_oralis, Moraxella_catarrhalis, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Staphylococcus_haemolyticus, Enterococcus_avium, Streptococcus_downei, Streptococcus_pyogenes, Streptococcus_salivarius_CAG_79, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_thermophilus, Eubacterium_rectale, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_invisus, Dialister_micraerophilus, Negativicoccus_succinicivorans, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Candida_orthopsilosis, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_concisus, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Serratia_marcescens, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:31 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:31 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:31 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:31 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:31 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:31 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 1, Actinobaculum_sp_oral_taxon_183
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 2, Actinomyces_georgiae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 3, Actinomyces_graevenitzii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 4, Actinomyces_hongkongensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 5, Actinomyces_johnsonii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 6, Actinomyces_massiliensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 7, Actinomyces_meyeri
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 8, Actinomyces_naeslundii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 9, Actinomyces_odontolyticus
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 10, Actinomyces_oris
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 11, Actinomyces_sp_HMSC035G02
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 12, Actinomyces_sp_HPA0247
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 13, Actinomyces_sp_ICM47
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 14, Actinomyces_sp_S6_Spd3
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 15, Actinomyces_sp_oral_taxon_170
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 16, Actinomyces_sp_oral_taxon_180
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 17, Actinomyces_sp_oral_taxon_181
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 18, Actinomyces_sp_oral_taxon_414
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 19, Actinomyces_sp_oral_taxon_448
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 20, Actinomyces_sp_oral_taxon_897
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 21, Actinomyces_viscosus
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 22, Aeriscardovia_aeriphila
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 23, Alloscardovia_omnicolens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 24, Bifidobacterium_dentium
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 25, Bifidobacterium_longum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 26, Scardovia_wiggsiae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 27, Corynebacterium_atypicum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 28, Corynebacterium_durum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 29, Corynebacterium_matruchotii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 30, Rothia_aeria
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 31, Rothia_dentocariosa
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 32, Rothia_mucilaginosa
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 33, Cutibacterium_acnes
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 34, Propionibacterium_acidifaciens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 35, Pseudopropionibacterium_propionicum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 36, Atopobium_deltae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 37, Atopobium_parvulum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 38, Atopobium_rimae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 39, Olsenella_profusa
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 40, Olsenella_scatoligenes
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 41, Olsenella_uli
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 42, Collinsella_aerofaciens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 43, Collinsella_intestinalis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 44, Collinsella_stercoris
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 45, Enorma_massiliensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 46, X.Collinsella._massiliensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 47, Adlercreutzia_equolifaciens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 48, Asaccharobacter_celatus
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 49, Denitrobacterium_detoxificans
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 50, Enterorhabdus_caecimuris
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 51, Slackia_exigua
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 52, Slackia_isoflavoniconvertens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 53, Rubrobacter_radiotolerans
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 54, Thermoleophilum_album
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 55, Porphyromonas_endodontalis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 56, Porphyromonas_somerae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 57, Prevotella_histicola
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 58, Prevotella_melaninogenica
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 59, Capnocytophaga_leadbetteri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 60, Kouleothrix_aurantiaca
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 61, Hydrogenibacillus_schlegelii
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 62, Gemella_bergeri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 63, Gemella_haemolysans
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 64, Gemella_morbillorum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 65, Gemella_sanguinis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 66, Brochothrix_campestris
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 67, Staphylococcus_argenteus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 68, Staphylococcus_aureus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 69, Staphylococcus_epidermidis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 70, Staphylococcus_schweitzeri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 71, Abiotrophia_defectiva
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 72, Abiotrophia_sp_HMSC24B09
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 73, Granulicatella_adiacens
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 74, Granulicatella_elegans
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 75, Enterococcus_faecalis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 76, Lactobacillus_rhamnosus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 77, Streptococcus_anginosus_group
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 78, Streptococcus_australis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 79, Streptococcus_cristatus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 80, Streptococcus_gordonii
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 81, Streptococcus_infantis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 82, Streptococcus_milleri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 83, Streptococcus_mitis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 84, Streptococcus_oralis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 85, Streptococcus_parasanguinis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 86, Streptococcus_peroris
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 87, Streptococcus_pneumoniae
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 88, Streptococcus_pseudopneumoniae
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 89, Streptococcus_salivarius
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 90, Streptococcus_sanguinis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 91, Streptococcus_sp_A12
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 92, Streptococcus_sp_F0442
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 93, Streptococcus_sp_HMSC034E03
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 94, Streptococcus_sp_HMSC067H01
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 95, Streptococcus_sp_HMSC071D03
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 96, Streptococcus_sp_HPH0090
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 97, Streptococcus_sp_M334
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 98, Streptococcus_vestibularis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 99, Eubacterium_brachy
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 100, Eubacterium_infirmum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 101, Eubacterium_nodatum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 102, Eubacterium_sulci
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 103, Mogibacterium_diversum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 104, Mogibacterium_pumilum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 105, Mogibacterium_timidum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 106, Lachnoanaerobaculum_saburreum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 107, Lachnospiraceae_bacterium_oral_taxon_096
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 108, Oribacterium_parvum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 109, Oribacterium_sinus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 110, Oribacterium_sp_oral_taxon_078
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 111, Shuttleworthia_satelles
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 112, Stomatobaculum_longum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 113, Peptostreptococcus_stomatis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 114, Bulleidia_extructa
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 115, Solobacterium_moorei
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 116, Limnochorda_pilosa
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 117, Veillonella_atypica
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 118, Veillonella_dispar
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 119, Veillonella_infantium
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 120, Veillonella_parvula
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 121, Veillonella_sp_T11011_6
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 122, Parvimonas_micra
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 123, Parvimonas_sp_oral_taxon_110
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 124, Parvimonas_sp_oral_taxon_393
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 125, Peptoniphilus_lacrimalis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 126, Fusobacterium_nucleatum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 127, Leptotrichia_wadei
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 128, Gemmata_obscuriglobus
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 129, Paludisphaera_borealis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 130, Achromobacter_denitrificans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 131, Achromobacter_ruhlandii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 132, Achromobacter_xylosoxidans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 133, Cupriavidus_sp
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 134, Sutterella_parvirubra
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 135, Anaerobiospirillum_thomasii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 136, Cardiobacterium_hominis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 137, Cardiobacterium_valvarum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 138, Alkalilimnicola_ehrlichii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 139, Escherichia_coli
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 140, Thiohalorhabdus_denitrificans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 141, Pseudomonas_aeruginosa_group
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 142, Fretibacterium_fastidiosum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 143, Acholeplasma_oculi
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 144, Candida_parapsilosis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 145, Bifidobacterium_breve
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 146, Nakamurella_silvestris
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 147, Porphyromonas_catoniae
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 148, Prevotella_jejuni
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 149, Prevotella_pallens
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 150, Prevotella_salivae
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 151, Tannerella_forsythia
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 152, Tannerella_sp_oral_taxon_HOT_286
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 153, Capnocytophaga_gingivalis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 154, Litorilinea_aerophila
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 155, Lactobacillus_gasseri
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 156, Sharpea_azabuensis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 157, Streptococcus_mutans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 158, Streptococcus_sp_oral_taxon_056
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 159, Oribacterium_asaccharolyticum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 160, Megasphaera_micronuciformis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 161, Neisseria_flavescens
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 162, Neisseria_subflava
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 163, Stenotrophomonas_maltophilia
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 164, Stenotrophomonas_pavanii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 165, Stenotrophomonas_rhizophila
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 166, Streptococcus_sobrinus
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 167, Streptococcus_viridans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 168, Eubacterium_saphenum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 169, Filifactor_alocis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 170, Eggerthia_catenaformis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 171, Leptotrichia_sp_oral_taxon_215
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 172, Haemophilus_parainfluenzae
## 2023-03-22 00:45:35 INFO::Counting total values for each feature
## 2023-03-22 00:45:35 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:35 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:35 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:35 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:35 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:35 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:35 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:35 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
cat("Number of differentially abundant bugs by each metadata in sputum")
## Number of differentially abundant bugs by each metadata in sputum
fit_data_spt$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                 8                 5                 1                 8 
##           molysis            qiaamp 
##                 8                 7
cat("\nAdding control data to both all and stratified MaAslin will help identify the actual contaminant\n")
## 
## Adding control data to both all and stratified MaAslin will help identify the actual contaminant
#Making significance table for figure
        # Define a function to make species names italicized
species_italic <- function(data) {
  names <- gsub("_", " ", rownames(data))
  names <- gsub("[]]|[[]", "", names)
  names <- gsub(" sp", " sp.", names)
  names <- gsub(" sp.", "* sp.", names)
  names <- gsub(" group", "* group.", names)
  names <- ifelse(grepl("[*]", names), paste("*", names, sep = ""), paste("*", names, "*", sep = ""))
  rownames(data) <- names
  data
}

#volcano plot of MaAslin with all data

ggplot(maaslin_all$results, aes(y = -log10(qval), x = coef, col = metadata)) +
        theme_classic(base_family = "serif") +
        #labs(tag = "A") +
        ggtitle("MaAslin with treatment types")+
        geom_point(size = 2) +
        xlab("MaAslin coefficient") +
        ylab("-log<sub>10</sub>(*q*-value)") +
        geom_hline(yintercept = 1, col = "gray") +
        geom_vline(xintercept = 0, col = "gray") +
        geom_richtext(aes( 4, 8, label = "*q*-value = 0.1, fold-change = 0", vjust = -1, fontface = 1), col = "grey", size = 3, family = "serif") +
        theme(legend.position = "top", axis.title.y = ggtext::element_markdown()) +
        scale_color_manual(values = c("#4daf4a",  "#984ea3", "#f781bf", "#377eb8", "#ff7f00", "#ffff33", "#a65628")) +
        guides(col = guide_legend(title = "Fixed effects", title.position = "top", nrow = 3))

# Make a significance table for each figure (top 20 taxa)


fit_data_spt$results
make_sig_table <- function(data) {
  sig_data <- spread(data$results[order(data$results$qval), c("feature", "metadata", "qval")], metadata, qval)
  sig_data$min <- apply(sig_data, 1, FUN = min)
  sig_data <- sig_data[order(sig_data$min),] %>% .[1:20,]
  sig_data[["feature"]] <- ifelse(sig_data[["feature"]] == "X.Collinsella._massiliensis", "[Collinsella]_massiliensis", sig_data[["feature"]])
  sig_data_italic <- sig_data %>% rownames_to_column(var = "-") %>% column_to_rownames(var = "feature") %>% species_italic %>% select(-c("-", "min", "log10.Final_reads"))
  sig_data_sig <- ifelse(sig_data_italic < 0.1, "*", NA)
  return(list(data = sig_data, data_italic = sig_data_italic, data_sig = sig_data_sig))
}

fit_data_bal <- make_sig_table(fit_data_bal)
fit_data_ns <- make_sig_table(fit_data_ns)
fit_data_spt <- make_sig_table(fit_data_spt)

fit_data_ns$data
cat("\n\nStratified Masslin - significant of each but by treatment in BAL\n")
## 
## 
## Stratified Masslin - significant of each but by treatment in BAL
fit_data_bal$data_sig
##                                      benzonase host_zero lypma molysis qiaamp
## *Actinomyces radingae*                      NA        NA    NA      NA     NA
## *Bulleidia extructa*                        NA        NA    NA      NA     NA
## *Candida dubliniensis*                      NA        NA    NA      NA     NA
## *Enterococcus avium*                        NA        NA    NA      NA     NA
## *Finegoldia magna*                          NA        NA    NA      NA     NA
## *Gemella morbillorum*                       NA        NA    NA      NA     NA
## *Peptostreptococcus* sp. MV1                NA        NA    NA      NA     NA
## *Peptostreptococcus stomatis*               NA        NA    NA      NA     NA
## *Rothia mucilaginosa*                       NA        NA    NA      NA     NA
## *Staphylococcus epidermidis*                NA        NA    NA      NA     NA
## *Staphylococcus schweitzeri*                NA        NA    NA      NA     NA
## *Varibaculum cambriense*                    NA        NA    NA      NA     NA
## *Actinomyces naeslundii*                    NA        NA    NA      NA     NA
## *Corynebacterium tuberculostearicum*        NA        NA    NA      NA     NA
## *Gemella asaccharolytica*                   NA        NA    NA      NA     NA
## *Gemella haemolysans*                       NA        NA    NA      NA     NA
## *Gemella sanguinis*                         NA        NA    NA      NA     NA
## *Negativicoccus succinicivorans*            NA        NA    NA      NA     NA
## *Sutterella parvirubra*                     NA        NA    NA      NA     NA
## *Prevotella oris*                           NA        NA    NA      NA     NA
cat("\n\nStratified Masslin - significant of each but by treatment in nasal swab\n")
## 
## 
## Stratified Masslin - significant of each but by treatment in nasal swab
fit_data_ns$data_sig
##                                benzonase host_zero lypma molysis qiaamp
## *Cupriavidus* sp.              "*"       "*"       "*"   "*"     NA    
## *Streptococcus oralis*         NA        NA        NA    "*"     NA    
## *Staphylococcus epidermidis*   NA        NA        "*"   NA      NA    
## *Sutterella parvirubra*        NA        "*"       "*"   NA      NA    
## *Aeriscardovia aeriphila*      NA        NA        NA    NA      NA    
## *Rothia mucilaginosa*          NA        "*"       NA    "*"     NA    
## *Corynebacterium atypicum*     NA        NA        NA    NA      NA    
## *Propionibacterium namnetense* NA        NA        NA    NA      NA    
## *Slackia isoflavoniconvertens* NA        NA        NA    NA      NA    
## *Cutibacterium avidum*         NA        NA        NA    NA      NA    
## *Staphylococcus aureus*        NA        NA        NA    NA      NA    
## *Staphylococcus argenteus*     NA        "*"       NA    NA      NA    
## *Collinsella massiliensis*     NA        NA        NA    NA      "*"   
## *Limnochorda pilosa*           NA        NA        NA    "*"     NA    
## *Corynebacterium accolens*     NA        NA        "*"   NA      NA    
## *Malassezia restricta*         "*"       NA        NA    NA      NA    
## *Collinsella intestinalis*     NA        "*"       NA    NA      NA    
## *Finegoldia magna*             NA        "*"       NA    NA      NA    
## *Gemella haemolysans*          NA        NA        NA    NA      NA    
## *Thermoleophilum album*        NA        NA        NA    NA      NA
cat("\n\nStratified Masslin - significant of each but by treatment in Sputum\n")
## 
## 
## Stratified Masslin - significant of each but by treatment in Sputum
fit_data_spt$data_sig
##                                  benzonase host_zero lypma molysis qiaamp
## *Streptococcus infantis*         "*"       NA        "*"   NA      NA    
## *Gemella haemolysans*            "*"       NA        "*"   NA      NA    
## *Eubacterium sulci*              "*"       "*"       NA    "*"     "*"   
## *Granulicatella elegans*         "*"       NA        "*"   NA      NA    
## *Rothia dentocariosa*            "*"       NA        "*"   NA      NA    
## *Streptococcus* sp. F0442        "*"       NA        "*"   NA      NA    
## *Bifidobacterium longum*         "*"       "*"       NA    "*"     "*"   
## *Pseudomonas aeruginosa* group.  NA        NA        NA    "*"     "*"   
## *Neisseria flavescens*           "*"       NA        "*"   NA      NA    
## *Rothia mucilaginosa*            NA        NA        "*"   NA      NA    
## *Actinomyces naeslundii*         NA        "*"       NA    "*"     "*"   
## *Denitrobacterium detoxificans*  NA        "*"       NA    "*"     "*"   
## *Actinomyces* sp. oral taxon 897 NA        "*"       NA    "*"     "*"   
## *Prevotella melaninogenica*      NA        NA        NA    "*"     "*"   
## *Gemmata obscuriglobus*          NA        NA        NA    "*"     NA    
## *Streptococcus gordonii*         NA        NA        "*"   NA      NA    
## *Corynebacterium durum*          NA        NA        NA    NA      NA    
## *Actinomyces* sp. HPA0247        NA        NA        NA    NA      NA    
## *Streptococcus sanguinis*        NA        NA        NA    NA      NA    
## *Actinomyces massiliensis*       NA        NA        NA    NA      NA
=======

>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229

Results After adding control data, MaAslin needs to be reanalyzed. Adding controls (mock communities) for each treatment group will show more statistically valid results in y ~ log(final reads) + sample_type + treatment, (re = subject_id))

A6. Decontam - stratified by treatment

--> both stratified and nonstratified

input of DNA concentration: 16S qPCR data

https://github.com/benjjneb/decontam/issues/33

Ben Callahan: But in the more limited testing on qPCR data the method still seems to work, and other publications report strong patterns of inverse frequency of contaminants using qPCR data - which is the pattern the frequency method relies on.

Strategy:

2.4.1. run decontam for all samples (common contaminants, by extraction)

2.4.2. stratify decontam analysis per each treatment method (contaminants by depletion methods)

# Decontam package --------------------------------------------------------

# common contaminants across all the treatment methods
#Decontam - were there any contaminants?#

sample_data(phyloseq$phyloseq_rel)$is.neg <- grepl("Neg. ext.", sample_data(phyloseq$phyloseq_rel)$sample_type)
phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_rel, S.obs != 0)

#With all sampels
dec_f_all <- isContaminant(phyloseq_rel_nz, method="frequency", conc="DNA_bac_well")
dec_p_all <- isContaminant(phyloseq_rel_nz, method="prevalence", neg="is.neg", threshold=0.5)
dec_c_all <- isContaminant(phyloseq_rel_nz, method="combined", neg="is.neg", conc = "DNA_bac_well")

cat("decontam frequency - all sample")
## decontam frequency - all sample
dec_f_all %>% subset(.,.$contaminant)
cat("decontam prevalence - all sample")
## decontam prevalence - all sample
dec_p_all %>% subset(.,.$contaminant)
cat("decontam combined - all sample")
## decontam combined - all sample
dec_c_all %>% subset(.,.$contaminant)
#Stratified by sample type

cat("decontam prevalence - BAL")
## decontam prevalence - BAL
subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg. ext.")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)
cat("decontam prevalence - Nasal swab")
## decontam prevalence - Nasal swab
subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal swab", "Neg. ext.")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)
cat("decontam prevalence - Sputum")
## decontam prevalence - Sputum
subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg. ext.")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)
cat("decontam frequency - BAL")
## decontam frequency - BAL
subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg. ext.")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam frequency - Nasal swab")
## decontam frequency - Nasal swab
subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal swab", "Neg. ext.")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam frequency - Sputum")
## decontam frequency - Sputum
subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg. ext.")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam combined - BAL")
## decontam combined - BAL
subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg. ext.")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam combined - Nasal swab")
## decontam combined - Nasal swab
subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal swab", "Neg. ext.")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam combined - Sputum")
## decontam combined - Sputum
subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg. ext.")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)

A7. LM of function alpha diversity

--> both stratified and nonstratified
sample_data <- sample_data(phyloseq$phyloseq_path_rpkm) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))

cat("Species richness of all samples\n\n")
## Species richness of all samples
cat("S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -271.32533 63.94609 72.47916 -4.2430322 0.0000643
sample_typeNasal swab 66.71179 39.44578 32.64336 1.6912273 0.1003209
sample_typeSputum 142.80087 36.31835 41.80438 3.9319203 0.0003111
treatmentlyPMA 59.31768 31.18858 64.39962 1.9019036 0.0616583
treatmentBenzonase 96.59262 31.65626 66.97650 3.0512965 0.0032643
treatmentHost zero 129.84121 32.25492 67.36194 4.0254694 0.0001466
treatmentMolysis 150.32258 32.65321 67.58822 4.6036081 0.0000189
treatmentQIAamp 86.27368 32.71224 67.61999 2.6373516 0.0103578
log10(Final_reads) 52.77944 11.01007 70.54314 4.7937427 0.0000088
sample_typeNasal swab:treatmentlyPMA -23.81209 39.44179 64.87664 -0.6037273 0.5481301
sample_typeSputum:treatmentlyPMA -5.02253 39.40167 62.83470 -0.1274700 0.8989755
sample_typeNasal swab:treatmentBenzonase -80.25929 38.01990 65.40988 -2.1109811 0.0386006
sample_typeSputum:treatmentBenzonase -52.43718 38.65224 63.57018 -1.3566400 0.1796940
sample_typeNasal swab:treatmentHost zero -106.90271 36.81521 64.23327 -2.9037647 0.0050492
sample_typeSputum:treatmentHost zero -76.35866 38.92443 62.75019 -1.9617158 0.0542352
sample_typeNasal swab:treatmentMolysis -106.97855 38.54815 66.00383 -2.7751931 0.0071693
sample_typeSputum:treatmentMolysis -108.15918 39.34418 62.63915 -2.7490513 0.0077999
sample_typeNasal swab:treatmentQIAamp -75.62038 36.68939 64.06591 -2.0610968 0.0433576
sample_typeSputum:treatmentQIAamp -40.59820 38.59881 63.08680 -1.0517993 0.2969039 
#Species richness  - stratified 
cat("Species richness at NS\n\n")
## Species richness at NS
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal swab")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 68.35891 68.76371 28 0.9941131 0.3286817
treatmentlyPMA 13.88909 17.11543 28 0.8114955 0.4239262
treatmentBenzonase 20.95703 16.24183 28 1.2903126 0.2074965
treatmentHost zero 56.43052 18.37320 28 3.0713495 0.0047048
treatmentMolysis 51.56131 16.32822 28 3.1578029 0.0037880
treatmentQIAamp 54.41631 19.71859 28 2.7596444 0.0100871
log10(Final_reads) 12.25116 10.45340 28 1.1719787 0.2510816 
cat("Species richness at BAL\n\n")
## Species richness at BAL
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -640.83129 99.14027 19.00295 -6.4638848 0.0000034
treatmentlyPMA 14.46191 33.76511 18.36963 0.4283092 0.6734060
treatmentBenzonase 16.49562 35.65062 19.98741 0.4627021 0.6485757
treatmentHost zero 39.81066 36.93641 19.98154 1.0778162 0.2939563
treatmentMolysis 54.18626 37.78713 19.92262 1.4339874 0.1670859
treatmentQIAamp -10.73947 37.91284 19.91125 -0.2832674 0.7798961
log10(Final_reads) 124.09731 17.94627 17.34430 6.9149360 0.0000022 
cat("Species richness at sputum\n\n")
## Species richness at sputum
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -127.58303 198.98655 21.38568 -0.6411641 0.5282282
treatmentlyPMA 54.38241 28.77622 19.88708 1.8898385 0.0734460
treatmentBenzonase 44.29211 36.29072 20.41549 1.2204804 0.2361930
treatmentHost zero 53.75292 61.00733 20.95720 0.8810895 0.3882627
treatmentMolysis 42.48496 71.13867 21.03201 0.5972132 0.5567405
treatmentQIAamp 45.90423 52.95334 20.86328 0.8668807 0.3958683
log10(Final_reads) 52.61786 33.96577 21.23460 1.5491436 0.1361234 
#Shannon
cat("Shannon index\n\n")
## Shannon index
cat("Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_shannon ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.6277442 0.1923582 72.90417 3.2634128 0.0016771
sample_typeNasal swab 0.3782276 0.0948640 66.24283 3.9870501 0.0001694
sample_typeSputum 0.3038993 0.0968430 65.80070 3.1380607 0.0025449
treatmentlyPMA 0.2080648 0.0912575 59.31888 2.2799751 0.0262196
treatmentBenzonase 0.1074536 0.0924809 65.10258 1.1619011 0.2495177
treatmentHost zero 0.2663135 0.0943090 66.22861 2.8238390 0.0062621
treatmentMolysis 0.3493594 0.0955289 66.88889 3.6571063 0.0005034
treatmentQIAamp 0.0974529 0.0957099 66.98134 1.0182112 0.3122414
log10(Final_reads) -0.0167122 0.0335618 71.90795 -0.4979546 0.6200356
sample_typeNasal swab:treatmentlyPMA -0.2199998 0.1167864 63.60638 -1.8837795 0.0641674
sample_typeSputum:treatmentlyPMA -0.0942524 0.1154524 55.89187 -0.8163741 0.4177521
sample_typeNasal swab:treatmentBenzonase -0.0695042 0.1120531 63.47768 -0.6202792 0.5372946
sample_typeSputum:treatmentBenzonase 0.0023288 0.1130661 57.07248 0.0205967 0.9836392
sample_typeNasal swab:treatmentHost zero -0.2682060 0.1085902 60.90331 -2.4698926 0.0163329
sample_typeSputum:treatmentHost zero -0.2089075 0.1142346 56.03885 -1.8287583 0.0727592
sample_typeNasal swab:treatmentMolysis -0.3849402 0.1134802 64.38350 -3.3921345 0.0011901
sample_typeSputum:treatmentMolysis -0.2757981 0.1156639 56.27671 -2.3844786 0.0204979
sample_typeNasal swab:treatmentQIAamp -0.1846564 0.1083181 60.79440 -1.7047607 0.0933449
sample_typeSputum:treatmentQIAamp -0.0636306 0.1130449 56.24157 -0.5628795 0.5757540 
#Simpson
cat("Inverse Simpson index\n\n")
## Inverse Simpson index
cat("InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_invsimpson ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.5338586 0.3135385 70.96294 8.0814913 0.0000000
sample_typeNasal swab 0.2842938 0.1520494 68.03030 1.8697465 0.0658228
sample_typeSputum 0.3782501 0.1558509 70.95435 2.4270003 0.0177657
treatmentlyPMA 0.3397324 0.1563898 61.65671 2.1723442 0.0336834
treatmentBenzonase 0.0472378 0.1563836 68.23317 0.3020637 0.7635223
treatmentHost zero 0.2329645 0.1590151 69.37014 1.4650462 0.1474261
treatmentMolysis 0.3567617 0.1607861 70.00663 2.2188598 0.0297400
treatmentQIAamp 0.0582910 0.1610497 70.09361 0.3619443 0.7184825
log10(Final_reads) -0.1644187 0.0540196 <<<<<<< HEAD 64.44675 ======= 64.44674 >>>>>>> a74da024aec8928028c246dbcfe35b96e9691229 -3.0436866 0.0033807
sample_typeNasal swab:treatmentlyPMA -0.4423346 0.1983466 65.97618 -2.2301099 0.0291475
sample_typeSputum:treatmentlyPMA -0.3759508 0.1993215 57.44567 -1.8861530 0.0643330
sample_typeNasal swab:treatmentBenzonase -0.0381545 0.1903215 66.04230 -0.2004738 0.8417261
sample_typeSputum:treatmentBenzonase -0.1934003 0.1947030 58.93845 -0.9933096 0.3246198
sample_typeNasal swab:treatmentHost zero -0.3256749 0.1855479 62.96653 -1.7552062 0.0840880
sample_typeSputum:treatmentHost zero -0.4518719 0.1971485 57.67798 -2.2920381 0.0255722
sample_typeNasal swab:treatmentMolysis -0.5649464 0.1923302 67.03548 -2.9373781 0.0045336
sample_typeSputum:treatmentMolysis -0.5317988 0.1995233 57.91638 -2.6653467 0.0099491
sample_typeNasal swab:treatmentQIAamp -0.1812611 0.1851750 62.57751 -0.9788638 0.3314165
sample_typeSputum:treatmentQIAamp -0.3137393 0.1950063 57.95487 -1.6088678 0.1130797 
#BPI
cat("Beger Parker index\n\n")
## Beger Parker index
cat("BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(dbp ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.3440391 0.1258762 70.13993 2.7331535 0.0079332
sample_typeNasal swab -0.0940641 0.0609175 68.34780 -1.5441228 0.1271764
sample_typeSputum -0.1495148 0.0626267 71.73071 -2.3873954 0.0196055
treatmentlyPMA -0.1360672 0.0637784 62.70891 -2.1334365 0.0368082
treatmentBenzonase -0.0137189 0.0634994 69.15451 -0.2160472 0.8295870
treatmentHost zero -0.0775264 0.0645073 70.21212 -1.2018236 0.2334698
treatmentMolysis -0.1514704 0.0651879 70.79134 -2.3235979 0.0230235
treatmentQIAamp -0.0147607 0.0652893 70.86956 -0.2260821 0.8217883
log10(Final_reads) 0.0701153 0.0215697 61.84125 3.2506433 0.0018652
sample_typeNasal swab:treatmentlyPMA 0.1871812 0.0806756 66.80034 2.3201722 0.0233951
sample_typeSputum:treatmentlyPMA 0.1784389 0.0814836 58.41893 2.1898742 0.0325394
sample_typeNasal swab:treatmentBenzonase 0.0181788 0.0774064 66.92125 0.2348489 0.8150431
sample_typeSputum:treatmentBenzonase 0.0966450 0.0795283 59.95203 1.2152280 0.2290433
sample_typeNasal swab:treatmentHost zero 0.1112523 0.0756141 63.85620 1.4713183 0.1461162
sample_typeSputum:treatmentHost zero 0.1766839 0.0805844 58.67215 2.1925319 0.0323198
sample_typeNasal swab:treatmentMolysis 0.2327702 0.0781693 67.88223 2.9777684 0.0040230
sample_typeSputum:treatmentMolysis 0.2309710 0.0815446 58.89704 2.8324495 0.0063138
sample_typeNasal swab:treatmentQIAamp 0.0548290 0.0754814 63.39830 0.7263905 0.4702729
sample_typeSputum:treatmentQIAamp 0.1291496 0.0796962 58.96349 1.6205246 0.1104560

A8. permanova of function beta diversity

--> both stratified and nonstratified
phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_path_rpkm, S.obs != 0 & sample_type %in% c("BAL", "Nasal swab", "Sputum"))

bray_perm_uni <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + treatment + subject_id,
                            data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)

bray_perm <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp + subject_id,
                            data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)

bray_perm_strata <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp,
                            data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F),
                            strata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F) %>% .$subject_id, permutations = 10000)

bray_perm_inter <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type * treatment + log10(Final_reads) + subject_id,
                                  data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

bray_perm_ns <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
                               data = subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab") %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

bray_perm_bal  <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "BAL"), method="bray") ~  lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
                                 data = subset_samples(phyloseq_rel_nz, sample_type == "BAL") %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

bray_perm_spt <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
                                data = subset_samples(phyloseq_rel_nz, sample_type == "Sputum") %>% sample_data %>% data.frame(check.names = F), permutations = 10000) 

cat("\nUnivariate analysis\n")
## 
## Univariate analysis
bray_perm_uni %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 1.557 0.344 53.603 0.00
log10(Final reads) 1 1.024 0.226 70.491 0.00
Treatment 5 0.124 0.027 1.703 0.14
Subject 10 0.766 0.169 5.276 0.00
Residual 73 1.060 0.234 NA NA
Total 91 4.530 1.000 NA NA 
cat("\nDetailed treatment\n")
## 
## Detailed treatment
bray_perm %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 1.557 0.344 53.603 0.000
log10(Final reads) 1 1.024 0.226 70.491 0.000
lyPMA 1 0.019 0.004 1.284 0.263
Benzonase 1 0.008 0.002 0.579 0.459
Host zero 1 0.008 0.002 0.564 0.460
Molysis 1 0.087 0.019 6.020 0.015
QIAamp 1 0.001 0.000 0.067 0.836
Subject 10 0.766 0.169 5.276 0.000
Residual 73 1.060 0.234 NA NA
Total 91 4.530 1.000 NA NA 
cat("\n Strata -detailed treatment\n")
## 
##  Strata -detailed treatment
bray_perm_strata %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 1.557 0.344 35.378 0.000
log10(Final reads) 1 1.024 0.226 46.525 0.000
lyPMA 1 0.019 0.004 0.847 0.424
Benzonase 1 0.008 0.002 0.382 0.521
Host zero 1 0.008 0.002 0.372 0.570
Molysis 1 0.087 0.019 3.973 0.054
QIAamp 1 0.001 0.000 0.044 0.964
Residual 83 1.826 0.403 NA NA
Total 91 4.530 1.000 NA NA 
cat("\n Interaction term \n")
## 
##  Interaction term
bray_perm_inter %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "sample_type:treatment" ~ 'Sample type X treatment',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
Sample type 2 1.557 0.344 66.198 0.000
Treatment 5 0.684 0.151 11.631 0.000
log10(Final reads) 1 0.463 0.102 39.412 0.000
Subject 10 0.766 0.169 6.515 0.000
Sample type X treatment 10 0.319 0.070 2.715 0.008
Residual 63 0.741 0.164 NA NA
Total 91 4.530 1.000 NA NA 
cat("\n Stratified - nasal swab \n")
## 
##  Stratified - nasal swab
bray_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
lyPMA 1 0.011 0.042 2.049 0.147
Benzonase 1 0.011 0.040 1.940 0.159
Host zero 1 0.009 0.033 1.615 0.205
Molysis 1 0.014 0.052 2.566 0.104
QIAamp 1 0.054 0.197 9.648 0.003
log10(Final reads) 1 0.028 0.103 5.051 0.027
Subject id 2 0.001 0.003 0.083 0.980
Residual 26 0.146 0.530 NA NA
Total 34 0.275 1.000 NA NA 
cat("\n Stratified - BAL \n")
## 
##  Stratified - BAL
bray_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
lyPMA 1 0.040 0.020 1.305 0.268
Benzonase 1 0.115 0.058 3.716 0.070
Host zero 1 0.027 0.014 0.886 0.370
Molysis 1 0.180 0.091 5.820 0.028
QIAamp 1 0.001 0.000 0.029 0.904
log10(Final reads) 1 0.603 0.306 19.535 0.001
Subject id 4 0.513 0.260 4.153 0.017
Residual 16 0.494 0.250 NA NA
Total 26 1.973 1.000 NA NA 
cat("\n Stratified - sputum \n")
## 
##  Stratified - sputum
bray_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")
Df SumOfSqs R2 F Pr(>F)
lyPMA 1 0.017 0.023 3.920 0.053
Benzonase 1 0.002 0.003 0.440 0.557
Host zero 1 0.064 0.089 15.072 0.001
Molysis 1 0.139 0.191 32.510 0.000
QIAamp 1 0.370 0.509 86.595 0.000
log10(Final reads) 1 0.032 0.045 7.588 0.011
Subject id 4 0.021 0.029 1.224 0.332
Residual 19 0.081 0.112 NA NA
Total 29 0.726 1.000 NA NA

A9. DA analysis for function, by sample type and treatment

--> both stratified and nonstratified
--> multiple levels?

No normalization will be conducted for RPKM data.

Thread: https://forum.biobakery.org/t/maaslin-with-shortbred-results-and-panphlan/3102/2

<<<<<<< HEAD

=======

>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229

Done.

#===============================================================================
#BTC.LineZero.Footer.1.1.0
#===============================================================================
#R markdown citation generator.
#===============================================================================
#RLB.Dependencies:
#   magrittr, pacman, stringr
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#BTC.Dependencies:
#   LineZero.Header
#===============================================================================
#Generates citations for each explicitly loaded library.
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
str_libraries <- c("r", str_libraries)
for (str_libraries in str_libraries) {
    str_libraries |>
        pacman::p_citation() |>
        print(bibtex = FALSE) |>
        capture.output() %>%
        .[-1:-3] %>% .[. != ""] |>
        stringr::str_squish() |>
        stringr::str_replace("_", "") |>
        cat()
    cat("\n")
}
## R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. We have invested a lot of time and effort in creating R, please cite it when using it for data analysis. See also 'citation("pkgname")' for citing R packages.
## Wickham H, Bryan J (2023). readxl: Read Excel Files_. R package version 1.4.2, <https://CRAN.R-project.org/package=readxl>.
## phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. Paul J. McMurdie and Susan Holmes (2013) PLoS ONE 8(4):e61217.
## Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
## Garbett SP, Stephens J, Simonov K, Xie Y, Dong Z, Wickham H, Horner J, reikoch, Beasley W, O'Connor B, Warnes GR, Quinn M, Kamvar ZN (2023). yaml: Methods to Convert R Data to YAML and Back_. R package version 2.3.7, <https://CRAN.R-project.org/package=yaml>. ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see 'help("citation")'.
## H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
## Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J (2022). vegan: Community Ecology Package. R package version 2.6-4, <https://CRAN.R-project.org/package=vegan>.
## Leo Lahti et al. microbiome R package. URL: http://microbiome.github.io
## Kassambara A (2023). ggpubr: 'ggplot2' Based Publication Ready Plots. R package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
## Simon Garnier, Noam Ross, Robert Rudis, Antônio P. Camargo, Marco Sciaini, and Cédric Scherer (2021). Rvision - Colorblind-Friendly Color Maps for R. R package version 0.6.2.
## Davis NM, Proctor D, Holmes SP, Relman DA, Callahan BJ (2017). "Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data." bioRxiv_, 221499. doi:10.1101/221499 <https://doi.org/10.1101/221499>.
## Auguie B (2017). gridExtra: Miscellaneous Functions for "Grid" Graphics. R package version 2.3, <https://CRAN.R-project.org/package=gridExtra>.
## Kassambara A (2023). ggpubr: 'ggplot2' Based Publication Ready Plots. R package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
## Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
## Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package: Tests in Linear Mixed Effects Models." Journal of Statistical Software, *82*(13), 1-26. doi:10.18637/jss.v082.i13 <https://doi.org/10.18637/jss.v082.i13>.
## Ooms J (2023). writexl: Export Data Frames to Excel 'xlsx' Format_. R package version 1.4.2, <https://CRAN.R-project.org/package=writexl>.
## Gonçalves da Silva A (2017). harrietr: Wrangle Phylogenetic Distance Matrices and Other Utilities. R package version 0.2.3, <https://CRAN.R-project.org/package=harrietr>.
## Mallick H et al. (2020). Multivariable Association in Population-scale Meta-omics Studies, http://huttenhower.sph.harvard.edu/maaslin2. To cite the MaAsLin 2 software, please use: Mallick H, Rahnavard A, McIver LJ (2020). MaAsLin 2: Multivariable Association in Population-scale Meta-omics Studies. R/Bioconductor package, http://huttenhower.sph.harvard.edu/maaslin2.
## Wilke C, Wiernik B (2022). ggtext: Improved Text Rendering Support for 'ggplot2'. R package version 0.1.2, <https://CRAN.R-project.org/package=ggtext>.
## Aphalo P (2022). ggpmisc: Miscellaneous Extensions to 'ggplot2'_. R package version 0.5.2, <https://CRAN.R-project.org/package=ggpmisc>.
## Auguie B (2017). gridExtra: Miscellaneous Functions for "Grid" Graphics. R package version 2.3, <https://CRAN.R-project.org/package=gridExtra>.
## Wood S, Scheipl F (2020). gamm4: Generalized Additive Mixed Models using 'mgcv' and 'lme4'. R package version 0.2-6, <https://CRAN.R-project.org/package=gamm4>. ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see 'help("citation")'.
## Hadley Wickham (2007). Reshaping Data with the reshape Package. Journal of Statistical Software, 21(12), 1-20. URL http://www.jstatsoft.org/v21/i12/.
## Zhu H (2021). kableExtra: Construct Complex Table with 'kable' and Pipe Syntax. R package version 1.3.4, <https://CRAN.R-project.org/package=kableExtra>.
## Yihui Xie (2023). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.42. Yihui Xie (2015) Dynamic Documents with R and knitr. 2nd edition. Chapman and Hall/CRC. ISBN 978-1498716963 Yihui Xie (2014) knitr: A Comprehensive Tool for Reproducible Research in R. In Victoria Stodden, Friedrich Leisch and Roger D. Peng, editors, Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595
#===============================================================================